tag:blogger.com,1999:blog-64120943544375305092024-03-28T08:53:21.680+09:00random dispersalRの統計ネタは自分自身の勉強を兼ね発信しています。作成者は統計の専門家ではないので自己責任でご参照ください。ご指摘も(応援も^^)歓迎します!Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.comBlogger40125tag:blogger.com,1999:blog-6412094354437530509.post-1985760942568665822024-01-02T11:56:00.011+09:002024-01-03T15:38:10.425+09:00情報量規準による多重比較のRコード<br /><br />AICなどの情報量規準を用いた多重比較の方法(Dayton 1998)を実装するRコードの紹介です。以前、共著論文で使用する機会があったのですが、久しぶりに掘り返して、だいぶスッキリしたコードへと改定できたので公開することにしました。手作業だと面倒くさかったのでだいぶ手間が省けるはずです。<br /><br />非常によく見かける解析方法に、ANOVAの後で多重比較(post-hoc tests)というコンボがあリます。いわゆる多重比較の方法は検定(頻度論)に基づいており、Type I errorの回避に特化しているのが気になっていました。一方で情報量規準に基づくモデル選択による解析を行うケースも多いですが、旧来の多重比較(頻度論)も併用するとダブルスタンダードになってしまいます。実際に査読で突っ込まれたこともあります。また多重比較は正規分布を仮定できるデータにしか適用できないことも制約になります。この点についてはHolmの方法など原理的に分布型に制約されない方法もありますが、組み合わせ数が多くなると検出力が低下する問題もあります。<br /><br />こういった問題をすっきりと解決できる、情報量規準に基づいたモデル選択による多重比較の方法があります。以前、北大・久保先生のスライドで知った手法ですが( https://kuboweb.github.io/-kubo/ce/2004/index.html )、下記のDayton(1998)が同様の手法(PCIC)を発表しており、根拠論文として参照可能です。特に手法名が提唱されているとは知らずに使っている例も多いかもしれませんが。また、名前が知られていないからこそ、統計ソフトにもRパッケージでも実装されてこなかったのかもしれません。<br /><br />Dayton CM (1998) Information criteria for the paired-comparisons problem. The American Statistician, 52(2):144-151. DOI: 10.1080/00031305.1998.10480554(ResearchGateからも取得可能)<br /><br />下記の10数行のコードで実装できました。lm, glm, glm.nbなどのモデル式を下記のpcic関数に食わせれば実行できるはずです。<br /><br /><span style="color: #04ff00;">require(partitions) </span><div># 組み合わせを表現する関数setpartsを使うのに必要なパッケージ、要インストール<br /><br /># pcic関数の定義、ここから:<br /><span style="color: #04ff00;">pcic <- function(Model, </span><span style="caret-color: rgb(4, 255, 0); color: #04ff00;">Data=Model$model, </span><span style="color: #04ff00;">Expl=names(Model$model)[2]) { </span></div><div><span style="color: #04ff00;"><span> </span></span><span style="color: #04ff00;"><span> </span>nLevels <- length(levels(Data[,Expl]))</span> # グループのリスト<span style="color: #04ff00;"><br /><span> </span>groupPTs <- apply(setparts(nLevels), 2, function(x) as.numeric(factor(x, levels=unique(x)))) <br /></span><span> </span># グルーピングの全パターン、setpartsの標準ではb,a,aなどが出るので正順に並べ替え<div><span style="color: #04ff00;"><span> </span>groupPTA <- apply(matrix(letters[1:nLevels][groupPTs], nrow=nLevels), 2, paste, collapse=" ") </span><div><span style="color: #04ff00;"><span> </span></span> # アルファベット変換&結合<span style="color: #04ff00;"><br /><span> </span>AICs <- c(AIC(update(Model, </span></div><div><span style="color: #04ff00;"><span> </span><span> </span>reformulate(c(setdiff(names(Data)[-1],Expl),1), </span><span style="color: #04ff00;">response = names(Data)[1]))), <span> </span></span> # 帰無モデル</div><div><span style="color: #04ff00;"><span> </span>sapply(2:ncol(groupPTs), function(i) { </span># 帰無モデル以外<span style="color: #04ff00;"><br /><span> </span>Data$groupLv <- factor(letters[1:nLevels][groupPTs[,i]][Data[,Expl]], levels=letters[1:nLevels]) </span></div><div><span style="color: #04ff00;"><span> </span></span># グループ分けの変数をDataに追加<span style="color: #04ff00;"><br /><span> </span>AIC(update(Model, reformulate(c(setdiff(names(Data)[-1],Expl),"groupLv"), </span></div><div><span style="color: #04ff00;"><span> </span><span> </span>response = names(Data)[1]), data=Data)) } ) ) </span># 各モデル<span style="color: #04ff00;"><br /><span> </span>comps <- data.frame(AIC=AICs, deltaAIC=AICs-min(AICs), Grouping=groupPTA)<br /><span> </span>print(paste(c(paste(levels(Data$Species), collapse=", "), "are labeled in this order.", <br /><span> </span><span> </span>Expl, "sharing the same alphabet do not differ in the", names(Data)[1]), collapse=" "))</span></div><div><span style="color: #04ff00;"><span> </span>return(comps[order(comps$AIC),]) }</span><br /># 関数定義ここまで(これを一度R上で実行するとカスタム関数として利用可能になる)</div><div><br />・lm, glm, glm.nbなどのモデルに適用可能です。それ以外でもModel$modelを実行してモデルに使用しているデータ(y, x)が返ってくるならば使用可能、offset項にも対応。</div><div>・注意点として、グルーピングのための説明変数は "~"のすぐ後ろの変数と見なす。異なる時は文字で指定する(この例の場合は"Species")<br /><div>・Data: Model$model # モデルに使用しているデータ:y, x。これで大丈夫なはずだが、自前で指定も可能</div></div><div><br /><br />実行例:<br /><span style="color: #04ff00;">data(iris)</span> # Sepal.Lengthが種によって異なるかについて多重比較を実行<span style="color: #04ff00;"><br />Model <- glm(Sepal.Length ~ Species, data=iris) <br /></span><span style="caret-color: rgb(4, 255, 0); color: #04ff00;">require(partitions)</span></div><div><span style="color: #04ff00;">pcic(Model)</span><br /></div></div><div><span style="color: #04ff00;"><br /></span></div><div><div>[1] "setosa, versicolor, virginica are labeled in this order. Species sharing the same alphabet do not differ in the Sepal.Length"</div><div> AIC deltaAIC Grouping</div><div>5 231.4520 0.00000 a b c</div><div>4 265.6359 34.18393 a b b</div><div>3 295.6778 64.22581 a a b</div><div>1 372.0795 140.62758 a a a</div><div>2 373.1310 141.67904 a b a</div></div><div><br /></div><div>一番上がベストモデルです。setosa: a, versicolor: b, virginica: c、のように異なるアルファベットで識別されており、種によってSepal.Lengthが異なることが示されました。このアルファベットによる識別は一般的な多重比較でよく見かける表現なので、多くの人が馴染みがあるかと期待します。</div></div><div><br /></div><div>結果を表形式で取り出すには、こんな感じでOKです。</div><div><span style="color: #04ff00;">write.csv(pcic(Model), "保存場所/ファイル名.csv")</span></div><div><br /></div><div>最後に、post-hoc testの用途について。実験系など個々のグループの意味が明確な場合に適した手法だと思いますが、その場合は頻度論の方が適していそうです。ここまで書いておいて何ですが、個人的にはpcicが必要なケースってそんなに多く無いと思っています。生態学系でよく見かける使用例の多くは回帰分析の方が適しているケースだし(説明変数が連続変数になっている等)、2way ANOVAなどの方が適しているケースも多いように感じます。何の仮定も置かずにpost-hoc testをやってみて、その結果を見てから演繹的に考察をするような使い方をしようとしていたら、少し立ち止まって考え直してみませんか。仮説を立ててそれを検証するための解析を行うという設計さえしていれば、より適切な解析方法が見つかるはずです。</div><div><br /></div><div><div>cf. lme4は今のところ非対応です。Data=Model@frameでデータを引っ張って来れそうですが、モデルのアップデートのところの改定がすぐにはできそうもなく。</div></div><div><br /></div>Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-62289139739441243282023-04-09T23:22:00.007+09:002023-04-25T11:20:58.474+09:00Rの常用速度を新旧macで比較:Apple Silicon M1 Max vs. Intel Core i9<p>今更ながら、macでのRの計算速度をIntel版とApple Siliconとで比較してみました。先行するテスト結果は検索するといくつか見つかるのですが、自分にとって普段よく使いそうな実用的用途で試してみたいと思いました。(2023.04.25追記:Intelも第12世代でCPU速度は飛躍的に向上して、以下の比較機種とは状況は異なっている模様)</p><p>比較をした計算環境は次の通りです。Rは登稿時の最新 ver. 4.2.3に揃えました。ハードウェアはintelとapple siliconとで実質1年半の差なので、ほぼCPUの差で比較できているかと思います。<br /><br />(1)macbook pro 14-2021(以下、MBP_M1mと略します)<br />Apple Silicon M1 Max 10コアCPU、32コアGPU 64GBメモリ<br /><br /><br />(2)macbook pro 16-inch 2019(以下、MBP_Intelと略します)<br />Intel Core i9 2.4GHz 8コアCPU 64GBメモリ<br /><br />(3)新旧macだけでもよかったのですが、Ubuntuを載せたカスタムPCも比較に加えました。<br />カスタムPC 2019年製(OS: Ubuntu 22.04.1 LTS)(以下、Intel_Ubuntと略します)<br />Intel Core i9-9980XE Extreme Edition 3.0-4.4GHz 18コアCPU 128GBメモリ<br /><br /></p><p>計算は次の4つのケースで実行し比較しました。概ねシングルコアの速度を試す形になっていますが、マルチコアを活用できるケースならば、単にコア数の多いマシンを使えばいいんじゃないかと思うので、制約の掛かりやすそうな事例をわざと選んでいます。</p><p><br /></p><b>### Case 1: Markov chain</b><div>いわゆるMCMC等の計算過程で動いているマルコフ連鎖、値の更新が必要なので一括処理できない計算の例の1つかと。</div><div><br /><span style="color: #04ff00;">x <- matrix(0, 10^6, 10) </span><br /># 10^6回分のマルコフ連鎖、10本の鎖を同時に計算<br /><span style="color: #04ff00;">system.time( for(r in 2:(nrow(x))) x[r,] <- 0.9*x[r-1,] + rnorm(10) )</span><br /><div><br /></div><div><br /></div><div><br /></div><div><div><b>### Case 2 raster_distance: raster画像の値あり全地点からの距離計算</b></div></div>都市や島からの距離の計算と思ったら良いです。<div><div><span style="color: #04ff00;">require(raster)</span></div><div><span style="color: #04ff00;">require(maptools)</span></div><div><span style="color: #04ff00;">data(wrld_simpl)</span></div><div><span style="color: #04ff00;">rgrid <- raster(ex=extent(c(-180,180,-90,90)), res=0.5)</span></div></div><div><div><span style="color: #04ff00;">rwrld <- raster::rasterize(wrld_simpl, rgrid)</span></div></div><div><div><span style="color: #04ff00;">system.time(rdst <- raster::distance(rwrld))</span></div></div><div>0.5度解像度の全世界の陸地からの距離を計算、かなり時間のかかるケース。</div><div>(ちなみにwrld_simplは気軽に呼び出しやすい粗い世界地図なのですが、{maptools}が2023年いっぱいで提供終了とのこと。高機能な{rnaturalearth}もよいですが...)</div><div><br /></div><div><br /></div><div><br /></div><div><div><b>### Case 3 raster_linear_model: raster stackへの回帰モデル適用</b></div></div><div>重ねたrasterの全値について回帰モデルを適用、calcのhelpにあるコードを高解像度にしただけのもの。単純な回帰モデル(lm)に限定(逆行列を利用した方法も紹介されていますが、応用しづらいので省きます)。</div><div><div><span style="color: #04ff00;">help(calc)</span></div><div><span style="color: #04ff00;">r <- raster(nrow=100, ncol=100) </span># ここを10倍にした</div><div><span style="color: #04ff00;">s1 <- lapply(1:12, function(i) setValues(r, rnorm(ncell(r), i, 3)))</span></div><div><span style="color: #04ff00;">s2 <- lapply(1:12, function(i) setValues(r, rnorm(ncell(r), i, 3)))</span></div><div><span style="color: #04ff00;">s1 <- stack(s1)</span></div><div><span style="color: #04ff00;">s2 <- stack(s2)</span></div><div><br /></div><div># regression of values in one brick (or stack) with another</div><div><span style="color: #04ff00;">s <- stack(s1, s2)</span></div><div># s1 and s2 have 12 layers; coefficients[2] is the slope</div><div><span style="color: #04ff00;">fun <- function(x) { lm(x[1:12] ~ x[13:24])$coefficients[2] }</span></div><div><span style="color: #04ff00;">system.time(x1 <- calc(s, fun))</span></div></div><div><br /></div><div><br /></div><div><br /></div><div><b>### Case 4: Stan (cmdstan)</b></div><div>コマンドライン版のStanで単純な状態空間モデルの計算、chain毎に異なるスレッドを割り当てているので、他の例と異なり、少しだけマルチコアの恩恵があるはず。コードは長いので、後に載せておきます。</div><div><br /></div><div><br /></div><div><br /></div><div>こちらが結果になります。左に行くほど高速です。</div><div><br /></div><div><div><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/a/AVvXsEgnQ3f7sn7OcAncvvKJqiI1eLEGJG0y1bjfwzf-wTLON2G544mThZJRRq7yq20cEJcA09nPoSl8uTrpzQ1NayOMVBGnn-3D83JOKlMuWW2j3IJd84A5bM6JHPOzokfuUJtv_uZy7E_p3zOHUVDiLf0TXUuPJPPkiejRAL9iCdQanK8nUUgwS1pYlfHO" style="margin-left: 1em; margin-right: 1em;"><img alt="" data-original-height="700" data-original-width="700" height="557" src="https://blogger.googleusercontent.com/img/a/AVvXsEgnQ3f7sn7OcAncvvKJqiI1eLEGJG0y1bjfwzf-wTLON2G544mThZJRRq7yq20cEJcA09nPoSl8uTrpzQ1NayOMVBGnn-3D83JOKlMuWW2j3IJd84A5bM6JHPOzokfuUJtv_uZy7E_p3zOHUVDiLf0TXUuPJPPkiejRAL9iCdQanK8nUUgwS1pYlfHO=w557-h557" width="557" /></a></div><div><br /></div><div><b><br /></b></div><div><b>概ねMBP_M1mはMBP_intelの2倍弱ほどの速度があることが分かりました。</b>計算時間の掛かるraster_distanceで2.33倍と最も差が大きいのは嬉しい結果です。マルチコアの恩恵があるStanでは1.5倍とさほどの差になっていないのも予想通りではあります。またOSの違いはあまり速度に関係ないようでした。</div><div><br /></div><div>なお、マルコフ連鎖については、2012年のマシンで5.457秒という手元の記録が残っているので(Late 2012 iMac Intel Core i7 3.4GHz クアッドコア 32GB)、倍速いというのがどれほどの進化かよく分かるかと思います。</div><div><br /></div><div>ちなみに、関連してCase 2, 3ではrasterパッケージの後継?terraパッケージでも比較をしてみたのですが、ほぼ変わりはなかったです。もしかすると連繋が進んでいたりして内部でterraを利用しているかもしれませんが、あくまで想像です。</div><div><br /></div><div><br /></div><div>もうM2が登場していますが、今回のようなRのシングルコア計算での比較ならば、おそらく速度は10%向上程度かなと見ています。</div><div><br /></div><div>その他、Apple SiliconはGPUコアや機械学習用のニューラルエンジンを積んでいるので、用途次第では次元の違う速度になる期待もあります。高速なGPUを使用可能なRパッケージってまだあまり聞かないようには思いますが。</div><div><br /></div><div>あとはSSDがかなり高速になっているそうで、物理メモリが不足した時にSSD領域を仮想メモリとして使用しても、かつてほどの速度低下はしないみたいで、もうさほどメモリに拘らなくてもよいのかもしれません。</div><div><br /></div><div><br /></div><div><br /></div><div><b># Case 4のcmdstanの計算コード:</b></div><div><span style="color: #04ff00;"><br /></span></div><span style="color: #04ff00;">require(cmdstanr)</span></div><div><span style="color: #04ff00;">stan_program <- "</span></div><div><span style="color: #04ff00;">data {</span></div><div><span style="color: #04ff00;"> int T;</span></div><div><span style="color: #04ff00;"> vector[T] Y;</span></div><div><span style="color: #04ff00;">}</span></div><div><span style="color: #04ff00;">parameters {</span></div><div><span style="color: #04ff00;"> real mu_Y0;</span></div><div><span style="color: #04ff00;"> vector[T] mu_Y;</span></div><div><span style="color: #04ff00;"> real<lower=0> sigma;</span></div><div><span style="color: #04ff00;"> real<lower=0> SIGMA;</span></div><div><span style="color: #04ff00;">}</span></div><div><span style="color: #04ff00;">model {</span></div><div><span style="color: #04ff00;"> // state model</span></div><div><span style="color: #04ff00;"> mu_Y[1] ~ normal(mu_Y0, sigma);</span></div><div><span style="color: #04ff00;"> for(t in 2:T) {</span></div><div><span style="color: #04ff00;"> mu_Y[t] ~ normal(mu_Y[t-1], sigma);</span></div><div><span style="color: #04ff00;"> }</span></div><div><span style="color: #04ff00;"><br /></span></div><div><span style="color: #04ff00;"> // observation model</span></div><div><span style="color: #04ff00;"> for(t in 1:T) {</span></div><div><span style="color: #04ff00;"> Y[t] ~ normal(mu_Y[t], SIGMA);</span></div><div><span style="color: #04ff00;"> }</span></div><div><span style="color: #04ff00;">}"</span></div><div><span style="color: #04ff00;"><br /></span></div><div><span style="color: #04ff00;">modfile <- write_stan_file(stan_program)</span></div><div><span style="color: #04ff00;"><br /></span></div><div><span style="color: #04ff00;">mod <- cmdstan_model(modfile) </span># 一時ファイルへ書き出し</div><div><br /></div><div><span style="color: #04ff00;">data(Nile)</span></div><div><span style="color: #04ff00;">data <- list(Y=Nile, T=length(Nile))</span></div><div><span style="color: #04ff00;"><br /></span></div><div><span style="color: #04ff00;">system.time(fit <- mod$sample(</span></div><div><span style="color: #04ff00;"> data = data,</span></div><div><span style="color: #04ff00;"> </span><span style="color: #04ff00;">chains = 4,</span></div><div><span style="color: #04ff00;"> parallel_chains = 4, </span># 並列化するchainの数</div><div><span style="color: #04ff00;"> iter_warmup = 5000, </span>#burn-in数</div><div><span style="color: #04ff00;"> iter_sampling = 20000 </span>#butn-in後のsampling数</div><div><span style="color: #04ff00;"> </span><span style="color: #04ff00;">)</span></div><div><span style="color: #04ff00;">)</span></div></div><div><br /></div></div>Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-15788103703759902352021-09-30T22:08:00.006+09:002021-10-02T11:37:04.162+09:00問題を起こすプロセスを定期的に強制終了する<p> macの動作を重くしたり不安定にさせたりするプロセス(アプリケーション含む)はなるべく終了させるか排除したくなるものです。純正のものではWindowServer, Falcond, mdworkerなどが重いですが、システムに必要そうなので我慢です。しかし、非純正のプロセスならいっそ落としてしまえば動作が改善する場合があるかもしれません。</p><p>このところmacが非常に不安定になっていて、Zoom会議で落ちる、Rで負荷の掛かる処理をすると落ちる、フォルダを移動しただけでもなど、処理内容が重い方が落ちやすいけれど軽くても構わず落ちる(カーネルパニック)現象がひどい時は1〜2日に一度くらいの頻度で起きるという困った状況に陥っていました。これではまったく使い物になりません。もしやと思い、心当たりの非純正アプリケーションを朝・夜にアクティビティモニタから強制終了させてみたところ、それらのトラブルは全て解消されました。残念ながら該当のアプリを削除するわけにはいかなさそうなのですが、この定期的なリフレッシュを自動で実行する仕組みを構築してみました。以下のような手順です。Automatorとカレンダー.appを使用します。</p><p><b>1. まず対象プロセスを強制終了させるアプリケーションを作成する。</b></p><p>1-1.「アプリケーション」から「Automator.app」を立ち上げる。</p><p></p><div class="separator" style="clear: both; text-align: center;"><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjqnmzOOkIpsXVi8KlT6HTwpDzxgjY0Dyr8iKN38c0Ev0kSmD1bgwCnx0-v38jM19FmljFVNZd0SajzOxA3nUeG7jHIzt2VT4uYQ_T0b0bc5x3VUm0isGNNADzih_TE2TFOWJIgip8oogo/s1510/automator%25E9%2581%25B8%25E6%258A%259E.jpg" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" data-original-height="298" data-original-width="1510" height="122" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjqnmzOOkIpsXVi8KlT6HTwpDzxgjY0Dyr8iKN38c0Ev0kSmD1bgwCnx0-v38jM19FmljFVNZd0SajzOxA3nUeG7jHIzt2VT4uYQ_T0b0bc5x3VUm0isGNNADzih_TE2TFOWJIgip8oogo/w618-h122/automator%25E9%2581%25B8%25E6%258A%259E.jpg" width="618" /></a></div><div class="separator" style="clear: both; text-align: left;"><br /></div><div class="separator" style="clear: both; text-align: left;">1-2. すると次のような新規書類が現れるはずなので「ワークフロー」を選択する。</div><div class="separator" style="clear: both; text-align: center;"><br /></div></div><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgQ33g1aRNMEjRSXnXbXJGgsdxn-FmKq2ofQICV6h2tYZ-y8moMOVzit__djk7YaBcLHyytv0oVgHkF10KydqtyMDLGPECp3rtYe5SiTiwK-n0FWpfH47ur0H7WS2OMXwd-Q03uvYD6yZs/s778/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588+2021-09-30+11.02.58.png" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" data-original-height="555" data-original-width="778" height="350" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgQ33g1aRNMEjRSXnXbXJGgsdxn-FmKq2ofQICV6h2tYZ-y8moMOVzit__djk7YaBcLHyytv0oVgHkF10KydqtyMDLGPECp3rtYe5SiTiwK-n0FWpfH47ur0H7WS2OMXwd-Q03uvYD6yZs/w491-h350/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588+2021-09-30+11.02.58.png" width="491" /></a><div class="separator" style="clear: both; text-align: left;"><br /></div><div class="separator" style="clear: both; text-align: left;">1-3. 「アクション」の中から「シェルスクリプトを実行」を選び、右の枠の中へとドラッグアンドドロップする。</div><div class="separator" style="clear: both; text-align: center;"><br /></div><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjDWeKkfduo0dPh0B1Cl9Gli2BuGnnam0je1NZ6_ekJ_h75_9Dr8frKrBKEdwXZ7zkm6TJLApPqtI_R6NBuEwjMOfqiP5VVbuItMqBPKAUDPCbJD7OqKFhKQYn9w_OwySHILeBLxQsi720/s940/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588+2021-09-30+11.31.08.png" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" data-original-height="415" data-original-width="940" height="256" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjDWeKkfduo0dPh0B1Cl9Gli2BuGnnam0je1NZ6_ekJ_h75_9Dr8frKrBKEdwXZ7zkm6TJLApPqtI_R6NBuEwjMOfqiP5VVbuItMqBPKAUDPCbJD7OqKFhKQYn9w_OwySHILeBLxQsi720/w581-h256/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588+2021-09-30+11.31.08.png" width="581" /></a></div><div class="separator" style="clear: both; text-align: center;"><br /></div><div style="text-align: left;">1-4. シェルスクリプトの入力は、「pkill プロセス名」とする。これはプロセス名のプロセス(アプリケーションである場合もある)を強制終了させるコマンドである。「ユーティリティ」の「アクティビティモニタ」にリストアップされるプロセス名で記す。プロセス名は先頭から一部の文字で構わない。また大文字と小文字を区別する必要がある。例えば、対象としたいプロセス名が「Hogehoge1」と「Hogehoge2」の場合、「pkill Hoge」で両方が対象となる。「pkill hoge」だと対象と見なされず何も起こらない。<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjQZTHzTM4Dwc7GhUqkCooobCsrpWeaPMod0Hj8Nqijp_YNSpqyAIw2-MuGvYhkF1E9gs-cHdIHazX6R9jY-LyUFhMrmtVA_Wq_8piBmVfNH9fON3N7lJadlgvnheppvnXSAmZ9XptMFmw/s989/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588+2021-09-30+11.42.57.png" style="clear: left; display: inline; margin-bottom: 1em; margin-right: 1em; text-align: center;"><img border="0" data-original-height="303" data-original-width="989" height="177" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjQZTHzTM4Dwc7GhUqkCooobCsrpWeaPMod0Hj8Nqijp_YNSpqyAIw2-MuGvYhkF1E9gs-cHdIHazX6R9jY-LyUFhMrmtVA_Wq_8piBmVfNH9fON3N7lJadlgvnheppvnXSAmZ9XptMFmw/w579-h177/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588+2021-09-30+11.42.57.png" width="579" /></a></div><div class="separator" style="clear: both; text-align: center;"> </div><div class="separator" style="clear: both; text-align: center;"><br /></div><div class="separator" style="clear: both; text-align: left;">1-5. 「ファイル」メニューから「アプリケーション」として保存する(ワークフローではなく)。保存場所はどこでも良さそう。一度、手動で立ち上げて、ちゃんと狙い通り動作(対象プロセスを強制終了)するか確かめておいた方が良いだろう。</div><div class="separator" style="clear: both; text-align: left;"><br /></div><div class="separator" style="clear: both; text-align: left;"><br /></div><div class="separator" style="clear: both; text-align: left;"><b>2. 作成したアプリケーションを定期的に実行させる仕組みを作る。</b></div><div class="separator" style="clear: both; text-align: left;">6. 「カレンダー.app」の予定を作成する。アプリを実行したい時間を設定する。さらに「通知」の仕方として「ファイル(イベントの開始時刻)を開きます」を選択、開始時刻リストの下にある「カスタム」を選択する。</div><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiQcqn608XfpNMtoVMqQY75Uil0UZ_7zp-ONBViL8JZ4Vpd-j60Pbmq40_tIAPiufiNWclsLlZkJAYooZbsXDQPGtItZHlJVWYesQ-b1jeGUY7c9SBnkROZDqpSpmNlJMYgYv7s9ntWqp8/s412/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588+2021-09-30+18.47.28.png" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" data-original-height="412" data-original-width="337" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiQcqn608XfpNMtoVMqQY75Uil0UZ_7zp-ONBViL8JZ4Vpd-j60Pbmq40_tIAPiufiNWclsLlZkJAYooZbsXDQPGtItZHlJVWYesQ-b1jeGUY7c9SBnkROZDqpSpmNlJMYgYv7s9ntWqp8/s320/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588+2021-09-30+18.47.28.png" width="262" /></a></div><div class="separator" style="clear: both; text-align: center;"><br /></div><div style="text-align: left;">7. 「ファイルを開く」対象がデフォルトでは「カレンダー」のところを、「その他...」を選択し、先ほど作成したアプリケーションを選択する。</div><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgR-9U6tAJHAnx7KbkdZNQKoc8Evrn-7SMt9RWQDbEG8qh0EjWUDrYMSsrXNLGSx_oRD3rKu3DQhu6Wlr5zvJkqWTBb_iH1ripNlzXmcNdTf8JvtLkXTP05TR34tpGntTErYWUAr488jgU/s328/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588+2021-09-30+18.47.54.png" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" data-original-height="328" data-original-width="326" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgR-9U6tAJHAnx7KbkdZNQKoc8Evrn-7SMt9RWQDbEG8qh0EjWUDrYMSsrXNLGSx_oRD3rKu3DQhu6Wlr5zvJkqWTBb_iH1ripNlzXmcNdTf8JvtLkXTP05TR34tpGntTErYWUAr488jgU/s320/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588+2021-09-30+18.47.54.png" width="318" /></a></div><div class="separator" style="clear: both; text-align: left;"><br /></div><div class="separator" style="clear: both; text-align: left;">8. 「繰り返し」を毎日にすれば、決まった時刻に実行される。私の場合、朝用と夜用を作成し、1日に2回リセットさせるように設定した。カレンダー.appは落としておいても動作するようです。</div><div class="separator" style="clear: both; text-align: left;"><br /></div><div class="separator" style="clear: both; text-align: left;">難点といえば「カレンダー」に毎日この定期リセットのための予定が書き込まれてしまうことです。ただし予定のグループを専用に設けて非表示にしておけば見えなくすることができます。また、通知機能を利用しているので、定時のリセットの度に通知されるのも気になると言えば気になりますが、ちゃんと機能していることを意識はできます。</div><br /><div class="separator" style="clear: both; text-align: center;"><br /></div><br /><div class="separator" style="clear: both; text-align: center;"><br /></div><br /><div class="separator" style="clear: both; text-align: center;"><br /></div></div><br /><br /><br /><p></p>Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-20389321331032248212020-12-20T19:27:00.007+09:002020-12-21T11:19:37.787+09:00tidyverse版: 複数生物種の個体数集計データを縦に伸ばすor縮める(集計データ← →データフレーム変換)<p> (<a href="http://nhkuma.blogspot.com/2012/12/or.html">以前の記事</a>をtidyverseで作り直してみました。まだ旧来のRコードの方が馴染みがあるので、tidyverseに慣れている方からすると不格好かもしれません。)</p><p>多数の種が含まれる生物群集データを解析する時、エクセルなどにデータを整理していると、各種をインデックスにするか要因のひとつとして整理するか、ケースによって変わってくるかと思います。とくに各種を応答変数とするか説明変数に使用するかあたりで必要になってくるかと。完全なデータフレーム型にしておいた方が統計解析はやりやすいけれど、膨大な行数となってエクセルの限界にあっさりと到達したりもするし…(cf. エクセルのvlookup関数)。</p><br />こんな時に、一部集計データとデータフレーム型とを自在に一発変換できたら効率がよいと思い、やり方をまとめてみました。正直、入力の手間が一番省けるのは部分的には集計されたデータだったりしますよね(下記のdata1のような形)。<br /><br /><div>データの整形にはdplyrパッケージの関数などを使用します。tidyverseパッケージとして関連パッケージも一括でインストールが可能です。</div><div><br /></div><span style="color: lime;">require(tidyverse)</span> # 要インストール、tidyverseのパッケージ群が一括で呼び出される。<div><br /><div># 例えば、こんな種毎に集計されているデータがある時<div># (tidyverse標準のtibble形式で作成)<br /><span style="color: lime;">data1 <- tibble(</span><br /><span style="color: lime;">year = c(rep("y08", 3), rep("y09", 3)), </span><br /><span style="color: lime;">site = c(1:3, 1:3), </span><br /><span style="color: lime;">depth = seq(1, 5, length=6), </span><br /><span style="color: lime;">sp1 = c(6:1), sp2 = c(1:6), sp3 = c(0:5))</span><br /><br /># 出力するとこんな感じです。<chr>は文字型、<int>は整数型、<dbl>は連続変数型を示す。<br />> data1<br /><div># A tibble: 6 x 6</div><div> year site depth sp1 sp2 sp3</div><div> <chr> <int> <dbl> <int> <int> <int></div><div>1 y08 1 1 6 1 0</div><div>2 y08 2 1.8 5 2 1</div><div>3 y08 3 2.6 4 3 2</div><div>4 y09 1 3.4 3 4 3</div><div>5 y09 2 4.2 2 5 4</div><div>6 y09 3 5 1 6 5</div><br /><br /><br />########################################################<br /><div><br /></div><span style="color: #3d85c6;"># data1をデータフレーム型に縦に伸ばす場合</span><br /><span style="color: lime;">data2 <- data1 %>% pivot_longer(c(-year, -site, -depth), names_to="species", values_to="abundance")</span><div># pivot_longerは集計された種名の表を縦長に伸ばす。種名を表す新たな種をspecies, 個体数をabundanceとする。</div><div># 形式を保持したい変数は、c(-year, -site, -depth)のように記す(c()は必要)。ひとつだけならば-yearのように書けば良いが、複数の場合は-c(year, site, depth)ではないことに注意。一般的なRの記法と異なっており、厄介な点。</div><div># "%>%"はパイプと呼ばれる。パイプの前のデータを後の処理に渡す意味。いちいちdata1$yearのように書かなくていい、with関数やattach関数と似ている。可読性は良くなっているかもだけれど、デバッグがしづらい気もする。</div><div><br />> data2 # 出力してみます<br /><div># A tibble: 18 x 5</div><div> year site depth species abundance</div><div> <chr> <int> <dbl> <chr> <int></div><div> 1 y08 1 1 sp1 6</div><div> 2 y08 1 1 sp2 1</div><div> 3 y08 1 1 sp3 0</div><div> 4 y08 2 1.8 sp1 5</div><div> 5 y08 2 1.8 sp2 2</div><div> 6 y08 2 1.8 sp3 1</div><div> 7 y08 3 2.6 sp1 4</div><div> 8 y08 3 2.6 sp2 3</div><div> 9 y08 3 2.6 sp3 2</div><div>10 y09 1 3.4 sp1 3</div><div>11 y09 1 3.4 sp2 4</div><div>12 y09 1 3.4 sp3 3</div><div>13 y09 2 4.2 sp1 2</div><div>14 y09 2 4.2 sp2 5</div><div>15 y09 2 4.2 sp3 4</div><div>16 y09 3 5 sp1 1</div><div>17 y09 3 5 sp2 6</div><div>18 y09 3 5 sp3 5</div><br /><br />########################################################<br /># もう一段階、伸ばしてみましょう。一個体あたり一行というデータセットへ。多項ロジットなど、カテゴリカルな解析に便利そうです。<br /><br /><span style="color: lime;">data3 <- select(data2, -abundance) %>% slice(unlist(map2(1:nrow(data2), data2$abundance, rep)))</span></div><div><div># data2のabundanceの数に応じてrepで行数を増やし、sliceで行を選択。map2はおおよそmapplyに相当、2変数を取るsapply。slice(unlist(map2(,,,,)))はもう少しスッキリ書けないものだろうか。</div><br />> data3 <br /><div># A tibble: 57 x 4</div><div> year site depth species</div><div> <chr> <int> <dbl> <chr> </div><div> 1 y08 1 1 sp1 </div><div> 2 y08 1 1 sp1 </div><div> 3 y08 1 1 sp1 </div><div> 4 y08 1 1 sp1 </div><div> 5 y08 1 1 sp1 </div><div> 6 y08 1 1 sp1 </div><div> 7 y08 1 1 sp2 </div><div> 8 y08 2 1.8 sp1 </div><div> 9 y08 2 1.8 sp1 </div><div>10 y08 2 1.8 sp1 </div><div># … with 47 more rows</div><div><br /></div><br /><br /><br />########################################################<br /><div><br /></div># data3を<span style="color: #3d85c6;">一つ前の状態に戻してみます(観測地点あたりで集計</span>、といったところ)<br /><span style="color: lime;">data4 <- data3 %>% count(year, site, depth, species, name="abundance") </span></div><div># year, site, depth, species毎にデータ頻度をカウントし、カウント結果をabundance変数として追加<br /><br /></div><div>> data4<br /><br /><div># A tibble: 17 x 5</div><div> year site depth species abundance</div><div> <chr> <int> <dbl> <chr> <int></div><div> 1 y08 1 1 sp1 6</div><div> 2 y08 1 1 sp2 1</div><div> 3 y08 2 1.8 sp1 5</div><div> 4 y08 2 1.8 sp2 2</div><div> 5 y08 2 1.8 sp3 1</div><div> 6 y08 3 2.6 sp1 4</div><div> 7 y08 3 2.6 sp2 3</div><div> 8 y08 3 2.6 sp3 2</div><div> 9 y09 1 3.4 sp1 3</div><div>10 y09 1 3.4 sp2 4</div><div>11 y09 1 3.4 sp3 3</div><div>12 y09 2 4.2 sp1 2</div><div>13 y09 2 4.2 sp2 5</div><div>14 y09 2 4.2 sp3 4</div><div>15 y09 3 5 sp1 1</div><div>16 y09 3 5 sp2 6</div><div>17 y09 3 5 sp3 5</div><div><br /></div># sp3 のゼロ個体データだけは消えてしまったが、当然と言えば当然<br /><br /><br /><br /><br />########################################################<br /><div><br /></div># data4の<span style="color: #3d85c6;">種数部分を集計表状に横に伸ばす</span>(最初の状態=data1の形にする場合)</div><div><br /><span style="color: lime;">data5 <- data4 %>% pivot_wider(names_from=species, values_from=abundance, values_fill=0)<br /> </span><br /># pivot_widerは縦長データから横長データへ変換する、見出しとなるspeciesと対応する値abundanceを展開する。生じるNAはvalues_fill=0で置き換えた。<br /><br />> data5 # 出力してみます。<br /><div># A tibble: 6 x 6</div><div> year site depth sp1 sp2 sp3</div><div> <chr> <int> <dbl> <int> <int> <int></div><div>1 y08 1 1 6 1 0</div><div>2 y08 2 1.8 5 2 1</div><div>3 y08 3 2.6 4 3 2</div><div>4 y09 1 3.4 3 4 3</div><div>5 y09 2 4.2 2 5 4</div><div>6 y09 3 5 1 6 5</div><br /><br /><br />########################################################<br /><div># 出現する種のリストにひどく変動がある場合など、1つのセルの中に種名を羅列したくなるケースもあると思います。そういうデータからの変換例も追加しておきます。</div><br /><span style="color: lime;">data6 <- tibble(site=paste0("s", 1:10), month=rep(c(1:2), each=5), </span><br /><span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>species=c("スズメ", "","スズメ、ヒヨドリ、シジュウカラ", "ムクドリ、スズメ", "ヒヨドリ、スズメ", "スズメ、キジバト", "スズメ、ムクドリ", "スズメ、ヒヨドリ", "", "ムクドリ"))</span><br /><br /># 例えばこんなかんじのデータセット。空欄さえあります…<br /><div><div># A tibble: 10 x 3</div><div> site month species </div><div> <chr> <int> <chr> </div><div> 1 s1 1 "スズメ" </div><div> 2 s2 1 "" </div><div> 3 s3 1 "スズメ、ヒヨドリ、シジュウカラ"</div><div> 4 s4 1 "ムクドリ、スズメ" </div><div> 5 s5 1 "ヒヨドリ、スズメ" </div><div> 6 s6 2 "スズメ、キジバト" </div><div> 7 s7 2 "スズメ、ムクドリ" </div><div> 8 s8 2 "スズメ、ヒヨドリ" </div><div> 9 s9 2 "" </div><div>10 s10 2 "ムクドリ" </div><div><br /></div># 注:このような日本語データのでの入力の場合、普通に読み込もうとすると文字化けしやすいので、readxlパッケージのread_excel関数などで読み込むのがトラブルが少なくていいと思っています。csvでは機種依存文字やmac/windows間での文字コード問題が混乱のもとなので。<br /><br /></div># こんなデータでも次のようにすれば、通常の集計データに変換できます。data1の形態へ変換してみましょう。<br /><span class="Apple-tab-span" style="white-space: pre;"> </span><br /><span style="color: lime;">sptabs <- data6 %>%</span> # data6をベースに集計作業<div><span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>pull(species) %>% </span># 種名部分をベクトルとして抽出</div><div><span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>map(str_split, pattern="、") %>%</span> # "、"で文字列を分割</div><div><span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>map_dfr(table) %>%</span> # 要素毎に集計テーブルを作成</div><div><span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>map_dfr(replace_na, 0) %>%</span> # naを0で置き換え</div><div><span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>select(-starts_with("..."))</span> # 出現無しを削除</div><div><span style="color: lime;">data7 <- select(data6, -species) %>%</span> # 元データの見出しを抽出</div><div><span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>mutate(sptabs)</span> # 種集計テーブルと結合</div><br />data7 # 出力してみます(出来上がり)<br /><div><br /><div># A tibble: 10 x 7</div><div> site month スズメ シジュウカラ ヒヨドリ ムクドリ キジバト</div><div> <chr> <int> <table> <table> <table> <table> <table> </div><div> 1 s1 1 1 0 0 0 0 </div><div> 2 s2 1 0 0 0 0 0 </div><div> 3 s3 1 1 1 1 0 0 </div><div> 4 s4 1 1 0 0 1 0 </div><div> 5 s5 1 1 0 1 0 0 </div><div> 6 s6 2 1 0 0 0 1 </div><div> 7 s7 2 1 0 0 1 0 </div><div> 8 s8 2 1 0 1 0 0 </div><div> 9 s9 2 0 0 0 0 0 </div><div>10 s10 2 0 0 0 1 0 </div></div></div></div><div><br /></div></div></div>Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-79784091499040695162020-09-03T15:27:00.016+09:002023-04-10T11:28:44.320+09:00Zoom / Teamsでスライドショーにならない/スライドが進まない問題の解決:macで外部ディスプレイ使用時?<p>コロナ禍のもと、オンラインツールを用いてプレゼンをする機会は多いと思います。<span style="color: #3d85c6;">Zoom / Teamsでスライドショーにならない/スライドが進まない</span>ことがあり、探しても解決策が見つからなかったため(書いてあっても単に私が理解していなかっただけかもですが)ここで紹介しておきます。<span style="color: #3d85c6;">いずれも外部ディスプレイの使用時</span>のトラブル例です。macでと書いていますが、Windowsで同様のことが起こるのかどうか単に未確認です。</p><p><br /></p><p><b><span style="color: #3d85c6;">Zoom、Teamsの両方に共通の解決方法が分かってきたので、その手順を先に書いておきます:</span></b></p><p>0. とっさの場合には、外部ディスプレイを引っこ抜いてmac本体単体で共有をやり直すのが手っ取り早いです。以下は、外部ディスプレイを保持したまま正常に共有でスライドショーにするための対処法です。</p><p>1. まずPowerPointでスライドショーを開始します。するとメイン画面は発表者ツール、外部ディスプレイが全画面のプレゼンになっているはずです。</p><p>2. PowerPointの発表者ツールの左上に、ツールのボタンが並んでいますが、「ディスプレイの入れ替え」をクリックします。</p><p></p><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhtlXB4l6Ozr2Mz3EgHxFhTCmssK-xRll2y69-hFnYKl0xUs1F48lTzq6NUuabqffXxTpICtS-KVtBnOsq-9KUHqgiKu-oWLxvbM-ku5SBR6DhZIYwmPaPVGhJfRqSbEzgK7QBrei2phn0/" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img alt="" data-original-height="125" data-original-width="428" height="93" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhtlXB4l6Ozr2Mz3EgHxFhTCmssK-xRll2y69-hFnYKl0xUs1F48lTzq6NUuabqffXxTpICtS-KVtBnOsq-9KUHqgiKu-oWLxvbM-ku5SBR6DhZIYwmPaPVGhJfRqSbEzgK7QBrei2phn0/" width="320" /></a></div>すると、メイン画面が全画面のプレゼン、外部ディスプレイが発表者ツールに入れ替わるはずです。<br />(2023追記:いっそここで「スライドショーの使用」を選び、両方のディスプレイにスライドショーが表示されるようにするのも良い解決方法かと思います。相手方を待たせるよりはいいでしょう。)<br /><p></p><p><br /></p><p><br /></p><p>3. この状態でZoom / Teamsへ移動し、通常通りにファイル共有ボタンをクリックします。この状態ではPowerPointのファイルを選択する形ではなく、ウィンドウ1または2の選択する形で全画面プレゼンになっているメイン画面を選択します。なお、PowerPointですでにプレゼンモードになっている状態でZoom / Teamsへ移動するには、command + タブ(このマークのボタン ->| )を押すと移動先のアプリを選択できるようになります。</p><p>4. これで相手側にも正常にZoom / Teams上でスライドショーが表示されているはずです。</p><p>5. プレゼンが終了したら、先に共有を終了してから、次にスライドショーを終了することをお勧めします。先にスライドショーを終了すると、デスクトップをさらけ出すことになります。</p><p>6. 一度この手順を踏んでおくと、その後、同じ環境では設定が保存されているようです。別のファイルを開こうとした際にも、初めから外部ディスプレイに発表者ツールが現れました。ただし、いざリアルのプロジェクターを用いたプレゼンに戻る際には注意が必要そうです。</p><p><br /></p><p>参考までに、自分自身ではどういう現象が起きたのかを書いておきます。</p><p><span style="color: #3d85c6;">●macのZoomファイル共有でPowerPointプレゼンのスライドショーが進まない場合</span></p><p>手元では共有したPowerPointをスライドショーに切り替え、手元ではスライドショーが進んでいるのに、相手側からは"最初のページから進んでいない"と言われました。後で実験してみると、相手側ではスライドショーモードにすら切り替わっていないことが分かりました。スライドショーを切って、編集モードの状態でスライドを進めると、相手側でも進みます。</p><div><br /></div><p><span style="color: #3d85c6;">●macのTeamsファイル共有でPowerPointプレゼン全画面できない場合</span> </p><p>手元では共有したPowerPointをスライドショーに切り替え、手元ではスライドショーが進んでいるのに、相手側からは"編集モードのままだ"と言われました。ディスプレイ設定をミラーリング、拡張のいずれもダメ、PowerPoint側を全画面/プレゼンツールへ切り替えなど試しましたが、自分の手元でしか全画面プレゼンにならずに困りました。</p>Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-49044842216511223352020-08-30T15:43:00.030+09:002021-12-31T14:15:58.491+09:00{cmdstanr}: Stan高速コンパイル、{rstan}代替としても有用<p>前回、最近<a href="https://nhkuma.blogspot.com/2020/04/macr.html">mac版のRパッケージ、とくにコンパイル必要パッケージのインストールが困難になっているという記事</a>を書いた。R3.6ではひどい状況だったが、R4.0.2の現在、CRANの手順に従えば、おおよそ問題はなくなっている。残る大きな問題は、rstanの動作不良だ。手元のR4.0.2では、rstanのインストール自体はできても、モデルのコンパイルができなくなっている。Macだけでなく、最近Windowsでもrstanの問題が出ている模様。</p><p>対策方法を調べていったところ、cmdstanrというコマンドラインベース(コンソール?macならTerminal)のRパッケージに辿り着いた。rstanのようにいちいちRcppでコンパイルしないので動作が速いし、またバージョンアップの際の依存環境との相性問題も生じにくいようだ。mac, linux, windowsのいずれでも使用可能なようだ。コマンドラインベースとは言うものの、Rから全ての作業を完結できることが分かったので、使い勝手も悪くない。</p><p><a href="https://mc-stan.org/cmdstanr/articles/cmdstanr.html">cmdstanrの公式ページ</a>が分かりやすいので、とくに説明も要らなさそうだが、実行サンプルを載せておく。</p><p># インストール: まず以下をR上で実行(公式ページの記述そのまま)</p><p><span style="color: #04ff00;">install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))</span></p><p># または、開発版のインストール(<strike>上記ではコンパイル失敗するモデルがあったが、そういう場合でも開発版では成功している2020.09.17追記 → </strike>2021.09.30現在、通常インストールの最新版で十分にうまく行くようになった模様)</p><p><span style="color: #04ff00;"></span></p><p> # おそらくもう不要→ #<span style="color: #04ff00;"> devtools::install_github("stan-dev/cmdstanr")</span></p><p># 次に、以下でcmdstanrを呼び出し、</p><p><span style="color: #04ff00;">require(cmdstanr)</span></p><p># R上で、下記のインストール関数を実行(coresの数値は、PCのコアをいくつ使用するかの指定、そこまで時間は掛からない)</p><p><span style="color: #04ff00;">install_cmdstan(cores = 2)</span></p><p><span style="color: #04ff00;">cmdstan_path()</span> # インストールされるディレクトリを示す場合</p><p># 2回目以降の使用で謎エラーが出てモデルがコンパイルできない場合、いったんcmdstanrを削除しインストールし直すことで使えている。いい対処法ではないだろうけれど、今のところ他にうまい手が見つからない。</p><p><br /></p><p>####################################</p><p># 計算実行:まず、モデル式を記述(公式サンプルではファイルからの読み込み)、</p><p><span style="color: #04ff00;">stan_program <- "</span></p><p><span style="color: #04ff00;">data {</span></p><p><span style="color: #04ff00;"> int<lower=0> N;</span></p><p><span style="color: #04ff00;"> int<lower=0,upper=1> y[N];</span></p><p><span style="color: #04ff00;">}</span></p><p><span style="color: #04ff00;">parameters {</span></p><p><span style="color: #04ff00;"> real<lower=0,upper=1> theta;</span></p><p><span style="color: #04ff00;">}</span></p><p><span style="color: #04ff00;">model {</span></p><p><span style="color: #04ff00;"> y ~ bernoulli(theta);</span></p><p><span style="color: #04ff00;">}</span></p><p><span style="color: #04ff00;">"</span></p><p><span style="color: #04ff00;"><br /></span></p><p><span style="color: #04ff00;">modfile <- write_stan_file(stan_program)</span></p><p><span style="color: #04ff00;">mod <- cmdstan_model(modfile)</span> # 一時ファイルへ書き出し、この方が運用しやすいと思う</p><p><br /></p><p># あとはrstanと同様の計算実行(公式の計算例、そのまま)</p><p><span style="color: #04ff00;">data_list <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))</span></p><p><span style="color: #04ff00;"><br /></span></p><p><span style="color: #04ff00;">fit <- mod$sample(</span></p><p><span style="color: #04ff00;"> data = data_list,</span></p><p><span style="color: #04ff00;"> seed = 123,</span> # 乱数の種</p><p><span style="color: #04ff00;"> chains = 4,</span></p><p><span style="color: #04ff00;"> parallel_chains = 4,</span> # 並列化するchainの数</p><p><span style="color: #04ff00;"> refresh = 500</span> # コンソール上に経過を表示するiterationの更新間隔</p><p><span style="color: #04ff00;">)</span></p># <a href="https://mc-stan.org/cmdstanr/reference/model-method-sample.html">samplingの詳細設定の解説も充実</a><div><b>## 指定可能な代表的なもの</b></div><div># iter_warmup: burn-in数</div><div># iter_sampling: butn-in後のsampling数</div><div># thin: sampleを保存する間隔</div><div># init: 初期値、init=function() list(theta=runif(1))のように</div><div><br /></div><div><br /></div><div># 計算結果はtibble形式になっている</div><div><span style="color: #04ff00;">fit$summary() </span># 計算結果の要約、tibble形式</div><div><span style="color: #04ff00;">fit$sampler_diagnostics() </span># 収束診断</div><div><span style="color: #04ff00;">fit$draws() </span># sample値(次元:iteration, chain, parameter)、</div><div># "draws_array"形式だが、<span style="color: #04ff00;">matrix(fit$draws()[,,1])</span>あるいは<span style="color: #04ff00;">matrix(fit$draws("theta"))</span>のようにparameter名で指定し、parameter毎に取り出し可能</div><div><br /></div><div><br /></div><div># しかし、以下によってrstanの形式に変換することも可能</div><div><span style="color: #04ff00;">stanfit <- rstan::read_stan_csv(fit$output_files())</span></div><div><span style="color: #04ff00;"><br /></span></div># rstanの形式にできた方が都合が良いだろうが、rstanのインストールが全くできない場合はtibbleの方がありがたい場合もあるだろう。<div><br /></div><div># 推定結果を保存する場合、以下のように.RDSで保存するのが推奨されている。</div><div><div><span style="color: #04ff00;">fit$save_object(file = "fit.RDS")</span></div><div># 再読み込みの場合</div><div><span style="color: #04ff00;">fit2 <- readRDS("fit.RDS")</span></div></div><div><span style="color: #04ff00;"><br /></span></div><div><br /><br />#(2021.12.31補足)個々のchainに1つのCPU threadを割り当てる並列化は上記のようにparallel_chainsの指定をするだけなので容易だが、chainあたり複数のthreadを割り当てるには、モデルコード自体にも工夫が要るようだ。GPUの利用なども可能なようだが、いずれにしてもモデル構造自体の制約が大きく、またN数が非常に大きいなどのビッグデータには効果が大きいが、生態学でよく用いるような複雑なプロセスモデルでは恩恵が少なそうである。<a href="https://www.blogger.com/u/1/#"></a><a href="https://mc-stan.org/docs/2_28/cmdstan-guide/parallelization.html">https://mc-stan.org/docs/2_28/cmdstan-guide/parallelization.html</a></div>Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-42030019813132431392020-04-03T15:57:00.014+09:002021-10-02T11:32:57.948+09:00(2021.10月追記:もう大丈夫そう)当面要注意、mac版Rのアップデート<div>(2021.10.01追記)R4.1.1へアップグレードした様子では2020年に発生していた下記の数多くの問題は既に解消されたように見えます。またApple Silicon用のRバイナリーパッケージは既にCRANからダウンロード可能になっているのも確認しました。</div><div><br /></div><div>(2020.08.29追記)クリーンインストール状態のR4.0.2をセットアップ試みました。CRANに、ソースからのパッケージコンパイルをする場合の方法が追記されていて(<a href="https://cran.r-project.org/bin/macosx/">https://cran.r-project.org/bin/macosx/</a>)、それ通りでrstan以外は適正にパッケージがインストールできるようになった模様。Xcode, command line tools, XQuartz, GNU Fortran 8.2をインストールせよとあります。注意点は、GNU Fortran 8.2はリンク先でfor Marvericsとなっていますが、Catalinaにインストールしても今のところ問題は起きていません(cf. <a href="https://thecoatlessprofessor.com/programming/cpp/r-compiler-tools-for-rcpp-on-macos/">https://thecoatlessprofessor.com/programming/cpp/r-compiler-tools-for-rcpp-on-macos/</a>)。</div><div>Catalinaでrstanがうまく動作しない問題は深刻なようで、範囲も広くて把握できていないので、とりあえず放置しています。<a href="http://nhkuma.blogspot.com/2020/08/cmdstanrstanrstan.html">より問題の少なそうなrstan代替のパッケージ紹介もしておきました。</a></div><div><br /></div><div>まとめると、上記のようにcranのガイドに従って手順を踏みさえすれば、rstan以外のインストール問題は解決しているようです。それでも以前のように、Rインストーラひとつで終わりとはいかなくなっているので、ハードルが上がってしまっており残念なことです。</div><div><br /></div><div><br /></div>*************(以下は、過去の参考履歴として残しておきます)****************************<br />(2020.04.14追記)CRANからRcpp1.0.4.6が通常版として配布され、ソースからのコンパイルが失敗する問題はある程度解決しているかもしれません。しかし依然、バイナリー版を見つけられない問題は残っていますし、ソースのコンパイルが通常の方法のみで行けない問題も残っています(tidyverseの必須パッケージの1つxml2など)。<br />
<br />現在(2020.03.31に確認)<span style="color: #3d85c6;"><b>mac版のRをアップグレード後に、パッケージのインストールがうまくいかなくなるトラブルが生じるかもしれません。どれくらいの範囲、環境まで影響するかわからないのですが、CatalinaだけでなくEl Capitanでも発生しています。</b></span>沢山の追加パッケージをインストールする必要のある方、とくにC++、Fortranコンパイル系のパッケージを利用している方は、しばらく様子を見た方がいいかもしれません。実態がわからず、復旧に相当な労力を要しました。通常ならばmac版Rのパッケージの再インストールなど、install.packages関数のリストのコピペ一発で済むはずでしたが、今回は日単位の時間を浪費する羽目に合いました。最近Linux版Rのセットアップをしたところでしたが、それよりもはるかに困難でした。<br />
<br />
生じた問題は大きく分けて2つあります。<br />
<br />
<b><span style="color: #3d85c6;">(1)install.packages関数がバイナリーパッケージを見つけられなかった。</span></b><br />
バイナリーよりもソースが新しいversionである場合に、ソースパッケージが選択されることはこれまでも時折ありましたが、今回はCRANのアドレスは正しいのに、必ずソースパッケージを探しに行ってしまう現象が生じています。試しにネットワーク環境を変えてみても状況は変わらずで、また一方で同じネットワーク環境でも未アップデートのR(R3.5.1@El Capitan)では問題なくバイナリーを取得できました。エラーメッセージは以下のような感じに出てきました。<br />
<br />
<span style="color: lime;">install.packages("パッケージ名") </span>を実行すると、以下のようなエラーメッセージが出て、バイナリーをスキップしてソースファイルのインストールが試みられました:<br />
<br />
<span style="color: orange;">警告: リポジトリー https://cran.ism.ac.jp/bin/macosx/el-capitan/contrib/3.6 に対する索引にアクセスできません :</span><br />
<span style="color: orange;"> URL 'https://cran.ism.ac.jp/bin/macosx/el-capitan/contrib/3.6/PACKAGES' を開けません </span><br />
<span style="color: orange;"> ソースパッケージ ‘</span>パッケージ名<span style="color: orange;">’ をインストール中です </span><br />
<span style="color: orange;"> URL 'https://cran.ism.ac.jp/src/contrib/</span>パッケージ名<span style="color: orange;">.tar.gz' を試しています </span><br />
(この後、ソースのコンパイルが始まる。つまりバイナリーの取得には失敗するけれど、ソースは取得できているので、CRAN自体には接続できている)<br />
<br />
<span style="color: #3d85c6;"><b>(2)大半のRパッケージでソースのコンパイルに失敗した。</b></span><br />
強制的にソースを取得させられるのに、それがことごとく失敗しました。展開されたコード内に、<span style="color: orange;">'uuid_t'; 'uid_t'; </span>といった大量のエラーが生じてコンパイルできない現象が生じました。<br />
<br />
これらの現象は、El Capitan、クリーンインストール後のCatalina両方で発生しました(つまり私は原因がわからずに、いい機会だからとOSのクリーンインストールまでしてしまった!)。初めは、El CapitanのR3.6.1を3.6.3にアップしたときにこの症状が発生したので、いったんRを完全消去後に3.6.1を入れ直しました。しかし症状は改善せず。コンパイルに失敗しているので開発環境をクリアにする必要を感じて、この機会にCatalinaにアップしたら、まさかの再発に苦しめられたという顛末です。<br />
<br />
<br />
<span style="color: #3d85c6;"><b>****** 解決策 ******************************************************</b></span><br />
(1)については、さしあたってCRANにある対象パッケージのバイナリーファイルのアドレスをコピーしてきて実行すれば、とりあえずインストールできました。多くのユーザにとって、これが最短の解決法だと思います。パッケージ名を入れ替えるだけなら簡単ですが、バージョンがわからないので、"cran パッケージ名"でググるなどして調べる必要があります(ターミナルで、R CMD INSTALLで行けないかと思ったけれど、これはソースしか取得できない模様)。私はこっちの解決法に気づくのが遅れて、(2)に取り組んでしまいましたが...。<br />
# dplyrの場合の例、URLを用いたバイナリーからのインストールコード:<br />
<span style="color: lime;">install.packages("https://cran.r-project.org/bin/macosx/el-capitan/contrib/3.6/dplyr_0.8.5.tgz”) </span><br />
<br />
(2)はより深刻です。心折れかけたところでu_riboさんに相談したところ、関連するトピックを見つけて教えてもらいました(<a href="https://github.com/RcppCore/Rcpp/issues/1060">https://github.com/RcppCore/Rcpp/issues/1060</a>)。<br />
これを読み解くと、問題として、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はテスト版のようで、インストールするなら覚悟が必要な代物です。<br />
こちらの問題の対処はRcpp 1.0.4.4だけでなく、複数の対処が必要です。重要な手順はいくつか判明していますが、効果があったかわからないものも含んでいます。それらを検証できる腕もないので、他の方の環境で、同じやり方で解決できるのかどうかは不明です。何の保証もなくお勧めできない最終手段だと思います。<br />
<br />
(a)Rのインストールからやり直した方がよさそうなので、ターミナルから以下のコードでRを消去しました。"存在しない"とエラーが返ってくるものもありました。<br />
<span style="color: magenta;">rm -rf /Applications/R.app</span><br />
<span style="color: magenta;">sudo rm -rf /Library/Frameworks/R.framework</span><br />
<span style="color: magenta;">sudo rm /usr/bin/{R,Rscript}</span><br />
<span style="color: magenta;">sudo rm /usr/local/bin/R</span><br />
<span style="color: magenta;">sudo rm ~/library/R/</span><br />
<br />
(b)Rのパッケージファイル(.pkg)のインストール(<a href="https://cran.ism.ac.jp/bin/macosx/">https://cran.ism.ac.jp/bin/macosx/</a>)、現時点の最新版のR3.6.3を選択しました。Catalina用だけは分けられているので注意です。<br />
なお、こちら(<a href="https://ryanhomer.github.io/posts/build-openmp-macos-catalina-complete">https://ryanhomer.github.io/posts/build-openmp-macos-catalina-complete</a>)のおすすめに従い、カスタムインストールでTcl/TkとTexinfoを外しました。<br />
<br />
(c)開発環境Command Line Toolsのインストール。ターミナルで以下のコードを実行、Gバイト単位のファイルを取得するので時間が掛かります。<br />
<span style="color: magenta;">sudo xcode-select --install</span><br />
<br />
(d)ターミナルで以下のディレクトリを開き、SDKのversionを確認します:<br />
<span style="color: magenta;">open /Library/Developer/CommandLineTools/SDKs/ </span><br />
# 私のCatalina環境の場合、MacOSX10.15.sdkが既にインストールされていました。より下のversionの場合、別の対処が必要かもしれません(<a href="https://github.com/RcppCore/Rcpp/issues/1060">https://github.com/RcppCore/Rcpp/issues/1060</a>)。<br />
<br />
(e)(入れていなかったら)homebrewのインストール。ターミナルで以下のコードを打つ<br />
<span style="color: magenta;">ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"</span><br />
<br />
(f)Tcl/TkとTexinfo、その他のインストール(homebrewユーザへは、bとは別途、こうした方がよいと勧められていました)<br />
<span style="color: magenta;">brew install tcl-tk texinfo</span><br />
<span style="color: magenta;">brew install llvm libomp</span><br />
<span style="color: magenta;">brew install gcc</span><br />
<br />
(g)~/.R/Makevarsファイルの書き換え。こちらの4.の項目を実行。必要だったかどうか不明(<a href="https://thecoatlessprofessor.com/programming/cpp/r-compiler-tools-for-rcpp-on-macos/">https://thecoatlessprofessor.com/programming/cpp/r-compiler-tools-for-rcpp-on-macos/</a>)<br />
<br />
h)clang7.0、gfortran6.1の.pkgファイルのインストール(これは随所で勧められていました)<br />
<a href="https://cran.r-project.org/bin/macosx/tools/clang-7.0.0.pkg">https://cran.r-project.org/bin/macosx/tools/clang-7.0.0.pkg</a><br />
<a href="https://cran.r-project.org/bin/macosx/tools/gfortran-6.1.pkg">https://cran.r-project.org/bin/macosx/tools/gfortran-6.1.pkg</a><br />
<br />
(i)Rcpp1.0.4.6のインストール <strike>Rcpp1.0.4.4のインストール</strike><br />
<strike>以下のRコードで入りましたが、まだテスト版とのことなので、入れるなら覚悟して入れてください</strike>: 問題のあったversionのRcppはアップデートされ、通常のインストール方法でRcpp1.0.4.6がインストールされるようになっています。<br />
<span style="color: lime;"><strike>install.packages("Rcpp", repos="https://rcppcore.github.io/drat")</strike></span><br />
<br />
(j)tidyverseの必須パッケージの1つxml2のコンパイルを試みる際に以下のようなエラーが頻発するので、これも実行。こちらを参考(<a href="https://github.com/r-lib/xml2/issues/223">https://github.com/r-lib/xml2/issues/223</a>)。<br />
ちなみに、ANTICONF ERRORはLinux版のRでのパッケージコンパイルでよく見かけますが、ここを見ると対処方法のヒントが得られやすいです。macでこれをやらなければならないかと思うと気が重いですが。<br />
<br />
# ターミナルで、以下を実行<br />
<span style="color: magenta;">brew install pkg-config</span><br />
# Rで、以下を実行<br />
<span style="color: lime;">install.packages("xml2", configure.vars = c("INCLUDE_DIR=/Library/Developer/CommandLineTools/SDKs/MacOSX10.15.sdk/usr/include/libxml2"))</span><br />
<br />
<span style="color: orange;">* installing *source* package ‘xml2’ ...</span><br />
<span style="color: orange;">** package ‘xml2’ successfully unpacked and MD5 sums checked</span><br />
<span style="color: orange;">Found pkg-config cflags and libs!</span><br />
<span style="color: orange;">Using PKG_CFLAGS=-I/usr/include/libxml2</span><br />
<span style="color: orange;">Using PKG_LIBS=-L/usr/lib -lxml2 -lz -lpthread -licucore -lm</span><br />
<span style="color: orange;">------------------------- ANTICONF ERROR ---------------------------</span><br />
<span style="color: orange;">Configuration failed because libxml-2.0 was not found. Try installing:</span><br />
<span style="color: orange;"> * deb: libxml2-dev (Debian, Ubuntu, etc)</span><br />
<span style="color: orange;"> * rpm: libxml2-devel (Fedora, CentOS, RHEL)</span><br />
<span style="color: orange;"> * csw: libxml2_dev (Solaris)</span><br />
<span style="color: orange;">If libxml-2.0 is already installed, check that 'pkg-config' is in your</span><br />
<span style="color: orange;">PATH and PKG_CONFIG_PATH contains a libxml-2.0.pc file. If pkg-config</span><br />
<span style="color: orange;">is unavailable you can set INCLUDE_DIR and LIB_DIR manually via:</span><br />
<span style="color: orange;">R CMD INSTALL --configure-vars='INCLUDE_DIR=... LIB_DIR=...'</span><br />
<span style="color: orange;">--------------------------------------------------------------------</span><br />
<br />
(k)さらに、頻発する <span style="color: orange;">"make: gfortran: No such file or directory"</span>というエラーに対処するため、~/.R/Makevarsにパスを追記をする。こちらを参考( <a href="https://github.com/merliseclyde/BAS/issues/1">https://github.com/merliseclyde/BAS/issues/1 </a>)。これによって大多数のエラーが解消されました。<br />
<span style="color: magenta;">open ~/.R/Makevars</span><br />
# これによって開かれたファイルに以下のパスを書き込んで保存。置き換えでなく追記です。<br />
F77 = /usr/local/gfortran/bin/gfortran<br />
FC = /usr/local/gfortran/bin/gfortran<br />
<br />
(l)ここまで95%くらいのパッケージのインストールに成功しましたが、ここまで来て、パッケージのURLから直接取得する方法にたどり着いたので、これ以上の対処はせず。<br />
<br />
<br />
(2)まったくひどい対処方法だと思うので、現状では個々のRパッケージをダウンロードしてきてインストールする(1)の方法が無難だと思います。Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-88954967126528180172018-08-22T19:47:00.003+09:002018-08-22T23:04:35.045+09:00生態系の熱帯化:藻場が狭まり、サンゴ群集が拡大する渾身の論文がPNAS(米国科学アカデミー紀要)出版になりました(<a href="http://www.nies.go.jp/whatsnew/20180820/20180820.html">プレスリリースも出ました</a>):<br />
<a href="http://www.pnas.org/content/early/2018/08/14/1716826115"><u>Naoki H. Kumagai</u>, 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. <i>Proceedings of the National Academy of Sciences of the United States of America</i>.</a><br />
(<a href="https://doi.org/10.1073/pnas.1716826115">https://doi.org/10.1073/pnas.1716826115</a>)<br />
<br />
研究内容や説明についてはプレスリリースに譲りますが、以下は書き切れなかった、もう少し個人的な部分や苦労についてメモしておきます。<br />
<br />
研究開始から約5年、初投稿から1年9ヶ月も掛かった大作です(Nは1ヶ月でエディターキック、Sの査読まで進んだものの手が届かず、PNASで2回の改訂ののち受理)。これまでに培った、移動分散、メタ個体群の考え方、データベースの構築、多くの観察データを集めて用いる統計モデリング技術、ベイズ推定、GISなど、全てを込めました。<br />
<br />
この論文では、主に1950〜2015年の、数多くの地道な観察記録(439文献、22,253調査記録)をしらみつぶしに探索し、観察年・位置情報を付加した上で、45種の生物種ごとの日本国内の分布変化を網羅しました。私はおもに海藻と魚類の分布変化記録の収集と整備を担当しました(サンゴも一部だけ)。興味のある方は<a href="http://www.pnas.org/content/pnas/suppl/2018/08/15/1716826115.DCSupplemental/pnas.1716826115.sapp.pdf">Supporting Information</a>のFig S2を見ていただけると45種の個々の分布変化が分かります。<br />
お陰で日本のあらゆる海岸線の地名には詳しくなり、市町村合併や海岸線の変化(埋め立て)の歴史も垣間見ました。<br />
<br />
この論文ではさらに、藻場を構成する海藻や、より南方から分布を拡げてくる造礁サンゴ、海藻を食害する魚類の分布変化速度を指標として、藻場とサンゴ群集の分布変化を解析しています。気候変動影響のもとで、海流の輸送といった大スケールの物理的要因および魚類による海藻の食害といった生物間相互作用が組み合わさることで、観察された分布変化をモデルによってうまく再現できることを示しました。これらの解析・モデルを最終的な形態にもっていく過程では、共著者の皆さまに大きく助けていただきました。この研究を通じて、私自身のモデリング技術、GIS技術とその理解も大きく向上したので、実によい勉強機会になったと思います。<br />
<br />
PNASなどのshort formatスタイルは限られた字数の中に、必要なことを凝縮しつつ、しかも読みやすくしなければならないという、実際やってみると非常に難しい作業でした。字数が短い=簡単と考える人も居るかもしれませんが、記述的な新発見の論文でもない限り、short formatの論文の執筆はフルペーパーを書くよりはるかに難しい(あるいは慣れていないと難しい)作業だと感じました。これについても論文を改訂する過程で共著者の皆さまには大きく助けていただきました。自分自身の文章力も向上できたと思います。<br />
<br />
この論文で述べたように、今後の日本の温帯藻場は温暖化の直接的影響だけでなく、魚類やウニなどの摂食圧の影響がいっそう深刻になると予想されます。またサンゴの分布拡大も気候の変化速度よりも遅いため、海藻藻場もサンゴ群集も何かしら人の手を加えて保全する手立てを考えていく必要があると思います。この論文がそういう議論の必要性を感じる1つの切っ掛けになればと思います。<br />
<br />
なお、実際の海中では同じ海(例えば1回のダイビングの中で共存するくらいの範囲)に藻場とサンゴ群集が居合わせることはよくあるけれど、いざこれを写真に収めようとすると適した場所はなかなかありません。査読2回目の改訂の際に偶然、この論文で取り上げた全生物要素がバランスよく生息するシンボル的な海(宇和島市津島町田之浜)を訪れることができました。論文のFig.1A中央、プレスリリースの写真1中央の写真は、南方性コンブ類、温帯性ホンダワラ類、南方性ホンダワラ類、南方性サンゴの全要素が凝縮された奇跡の一枚です。<br />
(参照:田之浜の調査でお世話になった、<a href="http://aquagate01.blog33.fc2.com/blog-entry-556.html">愛媛ダイビングセンターさんのブログ</a>でも取り上げていただきました)<br />
<div style="text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhYqQArgS6Q2H35GL5JaPknK3CQPgg7nL53UdVOyIqtrulIq7-ACr0K5EH8C5EyJdu5xHRTDuzp3Ys9f9z4AAeFoxky9MXutTzbfKMWmOr3PQl5ZdjzQidz9HL4t-bT0MW7HMm4DIFP3_Q/s1600/Photo1_1.JPG" imageanchor="1"><img border="0" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhYqQArgS6Q2H35GL5JaPknK3CQPgg7nL53UdVOyIqtrulIq7-ACr0K5EH8C5EyJdu5xHRTDuzp3Ys9f9z4AAeFoxky9MXutTzbfKMWmOr3PQl5ZdjzQidz9HL4t-bT0MW7HMm4DIFP3_Q/s320/Photo1_1.JPG" width="222" /></a></div>
Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-52461026404846656272018-05-29T18:48:00.001+09:002020-03-03T14:52:57.199+09:00Rいろは:グラフ作成ggplot2編###### ggplot2とは? #################<br />
# 標準仕様でもそれらしい科学的グラフを作ってくれるRパッケージ<br />
# (plot関数の場合、相当なアレンジをしないと使い物にならない)<br />
# 色分けや複合グラフ、重ね合わせなどに強い、応用の利きやすいメリットがある<br />
# さらにcowplotパッケージを適用して、より学術用に適したグラフ作成を紹介します<br />
<br />
# 項や足し算形式でグラフ要素を追加する形態になっている。例えば、散布図、棒グラフのようなグラフ形式や、回帰直線など。<br />
# 項の指定方法などが(plot関数と違って)英単語を元にしているものが多く連想しやすい<br />
<br />
# <a href="http://www.rstudio.com/wp-content/uploads/2015/12/ggplot2-cheatsheet-2.0.pdf">使い方の早見表も参考になります</a><br />
<br />
<span style="color: lime;">install.packages("ggplot2", dep=T)</span> # 未インストールの場合<br />
<span style="color: lime;">install.packages("cowplot", dep=T) </span># 未インストールの場合<br />
<span style="color: lime;">install.packages("reshape", dep=T)</span> # 未インストールの場合<br />
<span style="color: lime;">require(ggplot2)</span> # 呼び出し<br />
<br />
###### まず定番のdata"iris"を用いて、originalのggplot2で描いてみます #########<br />
<br />
<span style="color: lime;">data(iris)</span> # iris呼び出し<br />
<span style="color: lime;">head(iris)</span> # irisの構成をチェック<br />
#Sepal.Length Sepal.Width Petal.Length Petal.Width Species<br />
#1 5.1 3.5 1.4 0.2 setosa<br />
#2 4.9 3.0 1.4 0.2 setosa<br />
#3 4.7 3.2 1.3 0.2 setosa<br />
#4 4.6 3.1 1.5 0.2 setosa<br />
#5 5.0 3.6 1.4 0.2 setosa<br />
#6 5.4 3.9 1.7 0.4 setosa<br />
<br />
<br />
<br />
##### 散布図 ####################<br />
<span style="color: lime;">ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width, fill=Species)) + geom_point(pch=21, size=2)</span><br />
<br />
# ggplot(data=iris, # データ元の指定<br />
# aes(x=Sepal.Length, y=Sepal.Width, # x, yに用いる変数の指定(aes()内で変数を指定)<br />
# fill=Species)) # Speciesによって異なる色を使用<br />
# + geom_point( # + geom_xxx: xxxの部分でグラフ種類を指定<br />
# pch=21, # pch=21: 塗りと枠に異なる色を適用可能(21~25) # 全種に適用する場合、aes外に記述<br />
# size=2) # プロットサイズ(少し大きくしました)<br />
<br />
# 以下でも同様のグラフになります(fill=Speciesを後半に記述)<br />
<span style="color: lime;">ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width)) + geom_point(aes(fill=Species), pch=21, size=2)</span><br />
<br />
# このままでは背景がうるさく、学術用に不向き<br />
<span style="color: lime;">require(cowplot)</span> # cowplotパッケージを呼び出し、同じようにグラフを作成します<br />
<span style="color: lime;">theme_set(theme_cowplot())</span> # このコードを通します(すると以降、cowplotスタイルに変更される)<br />
<br />
<span style="color: lime;">ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width, fill=Species)) + geom_point(pch=21, size=2)</span><br />
# 背景が白地のスッキリしたものに変わりました。文字のサイズもちょうどいいです。<br />
<br />
<br />
<br />
##### 箱ひげ図 ####################<br />
<span style="color: lime;">ggplot(data=iris, aes(x=Species, y=Sepal.Width, fill=Species)) + geom_boxplot()</span><br />
<br />
<br />
<br />
##### 棒グラフ+SD ####################<br />
# サンプルデータの準備<br />
<span style="color: lime;">MEAN1 <- with(iris, tapply(Sepal.Width, Species, mean))</span> # 各種の平均("対象データ", "グループ分け", "関数名")<br />
<span style="color: lime;">SD1 <- with(iris, tapply(Sepal.Width, Species, sd)) </span># 各グループの標準偏差<br />
<span style="color: lime;">iris2 <- data.frame(Species=names(MEAN1), me=MEAN1, sd=SD1)</span><br />
# グラフの作成<br />
<span style="color: lime;">ggplot(data=iris2, aes(x=Species, y=me, fill=Species)) + geom_bar(stat="identity") + geom_linerange(aes(ymax=me + sd, ymin=me - sd))</span><br />
<br />
<br />
<br />
##### 2要因の棒グラフ + SD ####################<br />
# サンプルデータの準備<br />
<span style="color: lime;">iris3 <- rbind(</span> # 1.2倍したサンプルデータを追加<br />
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>data.frame(Species=names(MEAN1), me=MEAN1, sd=SD1, treat="C"),</span><br />
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>data.frame(Species=names(MEAN1), me=1.2*MEAN1, sd=1.2*SD1, treat="T")<span class="Apple-tab-span" style="white-space: pre;"> </span>)</span><br />
# グラフの作成<span class="Apple-tab-span" style="white-space: pre;"> </span><br />
<span style="color: lime;">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))</span><br />
# position, dodge: 接したグラフを作成するためのkeyword<br />
<br />
<br />
<br />
##### ヒストグラム ####################<br />
<span style="color: lime;">ggplot(iris, aes(x=Sepal.Width, fill=Species)) + geom_histogram()</span><br />
<br />
<br />
<br />
##### 時系列データのグラフ ####################<br />
<span style="color: lime;">data(airquality) </span># data"airquality"を使用<br />
<span style="color: lime;">head(airquality) </span># データの中身をチェック<br />
#Ozone Solar.R Wind Temp Month Day<br />
#1 41 190 7.4 67 5 1<br />
#2 36 118 8.0 72 5 2<br />
#3 12 149 12.6 74 5 3<br />
#4 18 313 11.5 62 5 4<br />
#5 NA NA 14.3 56 5 5<br />
#6 28 NA 14.9 66 5 6<br />
<br />
# 時間形式に変換<br />
<span style="color: lime;">airquality$Time <- </span><span style="color: lime;"><span style="caret-color: rgb(0, 255, 0);">as.POSIXct(</span></span><span style="color: lime;">strptime(with(airquality, paste(2018, Month, Day, sep="-")), "%Y-%m-%d"))</span> # (2020.03.03少しコードを修正)<br />
<span style="color: lime;">ggplot(airquality, aes(x=Time, y=Ozone)) + geom_point(pch=21, cex=2, fill="purple") + geom_line()</span><br />
<br />
<br />
<br />
##### 時系列データのグラフ(複数系列) ####################<br />
<span style="color: lime;">require(reshape)</span> # 要インストール<br />
<span style="color: lime;">airquality2 <- melt(</span> # データ構成を変換<br />
<span class="Apple-tab-span" style="white-space: pre;"> </span><span style="color: lime;">airquality[,c("Ozone","Temp","Time")], id.var="Time")</span><br />
<span style="color: lime;">head(airquality2)</span><br />
#Time variable value<br />
#1 2018-05-01 Ozone 41<br />
#2 2018-05-02 Ozone 36<br />
#3 2018-05-03 Ozone 12<br />
#4 2018-05-04 Ozone 18<br />
#5 2018-05-05 Ozone NA<br />
#6 2018-05-06 Ozone 28<br />
<br />
# 2パネル構成<br />
<span style="color: lime;">base <- ggplot(airquality2, aes(x=Time, y=value)) </span># ごちゃごちゃしてきたので分割します<br />
<span style="color: lime;">poi <- geom_point(pch=21, cex=2, aes(fill=variable)) </span><br />
<span style="color: lime;">lin <- geom_line(aes(colour=variable)) </span><br />
<span style="color: lime;">base + poi + lin + facet_wrap(~ variable, scale="free_y", ncol=1) </span># variableでパネルを分割<br />
<br />
<br />
<br />
<br />
###### 散布図+回帰直線 ##################<br />
<span style="color: lime;">summary(iris)</span> # 値の幅をチェック<br />
#Sepal.Length Sepal.Width Petal.Length Petal.Width Species <br />
#Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100 setosa :50 <br />
#1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300 versicolor:50 <br />
#Median :5.800 Median :3.000 Median :4.350 Median :1.300 virginica :50 <br />
#Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199 <br />
#3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800 <br />
#Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500 <br />
<br />
<span style="color: lime;">md1 <- lm(Sepal.Length ~ Sepal.Width, iris[iris$Species=="setosa",])</span> # 1種に限定<br />
<span style="color: lime;">coef(md1)</span><br />
#(Intercept) Sepal.Width<br />
#6.5262226 -0.2233611<br />
<br />
# 散布図 + 回帰直線<br />
<span style="color: lime;">x.interval <- seq(2, 5, 0.1) </span># x(Sepal.Width)の最小〜最大を0.1刻みで細かく分割<br />
<span style="color: lime;">est.Y <- predict(md1, newdata=data.frame(Sepal.Width=x.interval))</span><br />
<span style="color: lime;">fitted <- data.frame(Sepal.Width=x.interval, Sepal.Length=est.Y)</span><br />
# fitted <- data.frame(Sepal.Width=x.interval, Sepal.Length=coef(md1)[1] + coef(md1)[2]*x.interval) # 上と同じ<br />
<span style="color: lime;">base2 <- ggplot(iris[iris$Species=="setosa",], aes(x=Sepal.Width, y=Sepal.Length)) </span><br />
<span style="color: lime;">poi2 <- geom_point(fill="turquoise", pch=21, size=2) </span><br />
<span style="color: lime;">lin2 <- geom_line(data=fitted, colour="turquoise", lwd=1) </span># ここは"data="を明記する必要あり<br />
<span style="color: lime;">base2 + poi2 + lin2</span><br />
<br />
<br />
<br />
##### 散布図 + 回帰直線 + 信頼区間 ##################<br />
<span style="color: lime;">est.Y2 <- predict(md1, newdata=data.frame(Sepal.Width=x.interval), interval="confidence")</span><br />
<span style="color: lime;">head(est.Y2)</span><br />
#fit lwr upr<br />
#1 4.019981 3.753101 4.286860<br />
#2 4.089030 3.839589 4.338470<br />
#3 4.158079 3.925981 4.390177<br />
#4 4.227128 4.012251 4.442004<br />
#5 4.296177 4.098369 4.493984<br />
#6 4.365226 4.184291 4.546160<br />
<br />
<span style="color: lime;">fitted2 <- data.frame(Sepal.Width=x.interval, Sepal.Length=est.Y2[,1], upperCI=est.Y2[,3], lowerCI=est.Y2[,2])</span><br />
<span style="color: lime;">band <- geom_ribbon(data=fitted2, aes(ymin=lowerCI, ymax=upperCI), alpha=0.5, fill="turquoise", linetype="blank")</span> # alpha: 半透明(0~1)<br />
<span style="color: lime;">base2 + poi2 + band + lin2</span> # 重ねたい順で<br />
<br />
<br />
<br />
##### 散布図 + 回帰直線(複数グループ) ##################<br />
<span style="color: lime;">x.interval3 <- seq(1,5,0.1)</span><br />
<span style="color: lime;">new.dat <- data.frame( </span># 種数分、繰り返す<br />
<span class="Apple-tab-span" style="white-space: pre;"> </span><span style="color: lime;">Species=rep(unique(iris$Species),each=length(x.interval3)), </span><br />
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>Sepal.Width=rep(x.interval3, length(unique(iris$Species))))</span><br />
<span style="color: lime;">md3 <- lm(Sepal.Length ~ Sepal.Width*Species, iris)</span> # 全種を対象、Speciesによる違い(ANCOVA)<br />
<span style="color: lime;">est.Y3 <- predict(md3, newdata=new.dat)</span><br />
<span style="color: lime;">fitted3 <- data.frame(Species=new.dat$Species, Sepal.Width=new.dat$Sepal.Width, Sepal.Length=est.Y3) </span># Species列を追加<br />
<span style="color: lime;">base3 <- ggplot(iris, aes(x=Sepal.Width, y=Sepal.Length)) </span><br />
<span style="color: lime;">poi3 <- geom_point(aes(fill=Species), pch=21, size=2) </span><br />
<span style="color: lime;">lin3 <- geom_line(data=fitted3, aes(colour=Species), lwd=1)</span> # ここは"data="を明記する必要あり<br />
<span style="color: lime;">base3 + poi3 + lin3</span><br />
<br />
<br />
<br />
<br />
##### 複数タイプのグラフの組み合わせ(cowplot利用) ##################<br />
<span style="color: lime;">g1 <- ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width)) + geom_point(aes(fill=Species), pch=21, size=2) </span><br />
<span style="color: lime;">g2 <- ggplot(data=iris, aes(x=Species, y=Petal.Width)) + geom_boxplot(aes(fill=Species))</span> # 箱ひげ図<br />
<span style="color: lime;">+ theme(axis.text.x = element_text(angle=90, vjust=0.5))</span> # label方向<br />
<br />
<span style="color: lime;">g12 <- plot_grid(g1, g2, labels="auto", align="h") </span><br />
# labels="auto":小文字でラベリング(大文字は"AUTO")<br />
# 直接入力の例:labels=c("a","b")<br />
# align="h":図の端を揃える("h", "v", "hv")<br />
<br />
#### グラフの調整(例:凡例を統一、サイズ調整)<br />
<span style="color: lime;">g1n <- g1 + theme(legend.position="none")</span> # 一方の凡例消去<br />
<span style="color: lime;">bd <- panel_border(colour=1, size=1) </span># 外枠追加<br />
<span style="color: lime;">g1n2 <- plot_grid(g1n + bd, g2+ bd, labels="auto", </span><br />
<span style="color: lime;"> rel_widths=c(1, 1.4), </span># 横幅を1:1.4に変更<br />
<span style="color: lime;"> scale=0.9,</span> # パネル間の間隔を広げる<br />
<span style="color: lime;">vjust=0.2)</span> # ラベル位置(a, b)の調整<br />
<br />
<br />
<br />
#### ファイルの保存 ############################## <br />
<span style="color: lime;">save_plot("directry/filename.pdf", g1n2)</span> # cowplotの関数、お任せで適度なサイズになる。指定も可能(2020.03.03: ggsaveからsave_plotへ変更)<br />
<br />
#### 水平線・垂線 ##############################<br />
<span style="color: lime;">abl <- geom_hline(yintercept=0, linetype="dotted")</span> # 水平 ("dotted": 点線)<br />
<span style="color: lime;">abl2 <- geom_vline(xintercept=0, linetype="dotted")</span> # 垂直<br />
<br />
<br />
#### 軸とラベルの調整 ############################<br />
<span style="color: lime;">x.ax <- scale_x_continuous(breaks=c(1:5), limits=c(1, 5), expand=c(0.001, 0.001), name="Sepal width (mm)")</span><br />
# ex.<br />
<span style="color: lime;">base3 + poi3 + lin3 + x.ax</span><br />
<br />
<br />
#### 色の指定 ##############################<br />
# グレースケール<br />
<span style="color: lime;">set.fillG <- scale_fill_grey(start=0.5, end=1)</span> # グレー 〜 白の場合(塗りつぶし色)<br />
# "fill"の部分を"colour"にすると線の色の指定("shape": 形状、<br />
<br />
<span style="color: lime;">set.fillC <- scale_fill_manual(values=c("turquoise", "salmon", "blue"))</span> # 塗りつぶし色<br />
<span style="color: lime;">set.colC <- scale_colour_manual(values=c("turquoise", "salmon", "blue"))</span> # 線の色<br />
# ex.<br />
<span style="color: lime;">base3 + poi3 + lin3 + set.fillC + set.colC</span><br />
<br />
#### cf. 色コードの選択 ##########################<br />
<a href="http://colorbrewer2.org/#type=sequential&scheme=BuGn&n=3">"見やすい"色の組み合わせ例の参照サイト</a>(ColorBrewer)<br />
fill="#2ca25f" のようにしてカラーコードで色の指定が可能<br />
<br />
<br />
<br />
#### 90° 回転(グラフを横向きにする)#######################<br />
<span style="color: lime;">+ coord_flip()</span><br />
<br />
<br />
#### 凡例の調整 ############################## <br />
<span style="color: lime;">+ theme(legend.position="none")</span> # 凡例消し<br />
<span style="color: lime;">+ theme(legend.position="bottom")</span> # ラベルの右揃え(cf. vjust)<br />
<br />
# 長すぎるラベルを改行する( \n が使える)<br />
<span style="color: lime;">scale_x_continuous(name = "Length\n (cm)")</span><br />
<br />
# 上付き、ギリシャ文字など<br />
<span style="color: lime;">expression(m^2) </span><br />
<span style="color: lime;">expression(Delta)</span> # Δ<br />
<br />
<br />
<br />
#### cowplotの調整 ##############################<br />
<span style="color: lime;">+ panel_border(colour=1, size=1) </span># defaultの線・文字色がなぜかグレーなので黒に変更<br />
<span style="color: lime;">+ theme_cowplot(font_size = 12)</span> # 文字が大きく感じる場合はサイズ指定<br />
<br />
# 日本語の利用(ただしアウトライン化されてしまう模様 )<br />
<span style="color: lime;">ggsave(..., family="Japan1GothicBBB",...) </span>Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-82464418313144976122017-03-30T18:09:00.004+09:002019-09-17T13:32:07.852+09:00Rで機種依存文字の混じったデータファイルを読み込む方法Rへ読み込む際に非常に厄介な機種依存文字(=環境依存文字)、<a href="http://nhkuma.blogspot.jp/2012/12/r_5.html">データを作成する際には絶対に避けてほしい代物ですが</a>、丸に数字など、日本のあらゆる文書で愛されており、根絶は非常に困難です。<br />
<br />
見渡せる程度の小さなファイルならば手作業で特定し、検索・置換で除去するのですが、巨大なファイルの場合、いったいどんな文字が悪さをしているのか見つけることさえできない場合がありますね。<br />
<br />
今回、到底手作業で対処できないサイズのファイルに取り組んでいた際にいい方法を思いついたので挙げておきます(小ネタです)。<br />
<br />
1) まず、機種依存文字のファイルは.xlsxファイルとして保存しておきます(ここでは"データファイル名.xlsx"とします)。<br />
<br />
2) 通常、read.csv関数などでデータファイルを読み込むところを、<span style="color: lime;">openxlsx</span>パッケージの<span style="color: lime;">read.xlsx</span>関数を使用します(要インストール、<span style="color: lime;">install.packages(</span>"openxlsx"<span style="color: lime;">)</span> などで)。<br />
<br />
3) 以下で読み込むことができます、fileEncodingの指定も特に不要の様子。<br />
<span style="color: lime;">data <- read.xlsx(</span>"データファイル名.xlsx"<span style="color: lime;">)</span><br />
<span style="color: lime;"><br /></span>
cf. もちろんファイル出力にも対応しています(outputというデータフレームを出力する場合の例)<span style="color: lime;"><span style="color: lime;">write.xlsx(</span></span>output, "出力ファイル名.xlsx"<span style="color: lime;"><span style="color: lime;">)</span></span><br />
<span style="color: lime;"><br /></span>
Excelファイルを直接読むのはまだ抵抗はあるのですが、気にしなければ汎用性のある方法なので便利です。あとはこのパッケージがちゃんと存続してくれることを祈るばかりですが。なお、他にもExcelファイルを読めるパッケージは複数存在しているようですが、他のは基本的にJavaに依存しているようなので却下しました。<br />
<br />
(2019.09.17 追記)<br />
列名に漢字が混じる場合にフリガナが追加されるエラーが出ていて対処方法が分からず。代替関数として、<span style="color: lime;">readxl</span>パッケージの<span style="color: lime;">read_excel</span>関数の方が良いかもしれません。こちらならこのエラーは起こらず、またJava不使用です。tidyverseの一部のようなので、将来性もあるかもしれません。一方、読み込むとデータフレームではなくtibleという新形式になります。tible非対応のパッケージもまだ多いので、通常のデータフレームに変換するには、<span style="color: lime;">as.data.frame</span>(tible_no_data)としてやればよいです。<br />
<br />
cf. 機種依存文字が混じっているファイルをread.csvで読もうとすると、"In scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings, : 入力コネクション 'ファイル名.csv' に不正な入力がありました "のようなエラーメッセージが出る。<br />
Rじゃなくとも、メールの送受信でも文字化けの温床となります。データ解析に限らず、トラブルの原因にしかなりません。Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-18524725626824995422016-12-21T13:38:00.002+09:002016-12-21T13:40:18.579+09:00macOS10.12 (Sierra) でWinBUGSを動かす (Wine, R2WinBUGS使用)メインマシンをクリーンインストールする羽目に会いました。今度こそWinBUGSは卒業したくはあるのですが、依然StanやJAGSでは通らないコードも残っておりまだ手放せずにいます。早くStanが離散パラメータを直接扱えるようになってほしい…<br />
<br />
以下は、いったんMountain lion (OSX10.8)をクリーンインストール → Sierraにアップグレードしたマシンにインストールした場合の報告例です。El Capitan(OSX10.11)の時とほぼ同様です。<br />
<br />
まず、XQuartzとHomebrew(パッケージ化されていないアプリを容易にインストールするための補助ツール?)を入れ、Homebrewを用いてWine(非Windows OS上でWin専用アプリケーションを実行する環境)をインストール、Wineのディレクトリ内にWinBUGSをインストールするという流れです。<br />
<br />
cf. 以下の「ターミナル」の使い方:<br />
「アプリケーション」→「ユーティリティ」にある「ターミナル」を立ち上げる。<br />
<br />
コンピュータ名:~ ユーザ名$<br />
<br />
このドルマーク $ の後にコマンドを打っていく。<br />
なお、インストールに関わるところでパスワードを求められるが、その都度、自分のアカウントのパスワードを入れる。<br />
(以降、パスワードを入れる作業は説明を省略)<br />
<div>
<br /></div>
<div>
<br />
<br />
<br /></div>
<div>
<以下、作業手順(もしかするとSierraでは1〜4の行程は不要かもしれない></div>
0)実行環境<br />
・Mountain lion (OSX10.8)をクリーンインストール → SierraにアップグレードしたiMac<br />
<br />
<br />
1) 下準備の開始、/usr/local/フォルダを作る、ロックをいったん外して操作をするという動作をするのですが、その前にシステムの基本的なセキュリティを一旦外します。推奨されていない動作だということをお忘れなく。<br />
<br />
リカバリモード(⌘+R を押しながら起動)で起動し、ターミナルを立ち上げる<br />
<br />
<br />
2) ターミナルの $ マークの直後に、下記のコードを打ち込む(セキュリティを外す作業)これはコピペでOK、以降も同様。<br />
<br />
<span style="color: magenta;">csrutil disable</span><br />
<br />
<br />
3) 通常の再起動をする<br />
<br />
<br />
4) ターミナルに下記を打ち込む(改行されて見えているだろうが、改行無しで打ち込む)<br />
<br />
<span style="color: magenta;">sudo mkdir /usr/local && sudo chflags norestricted /usr/local && sudo chown $(whoami):admin /usr/local && sudo chown -R $(whoami):admin /usr/local</span><br />
<br />
ただし、今回の実行環境では、このディレクトリは既に存在すると言われた。なので、この危なっかしい工程はもはや省略できるかもしれない。<br />
<br />
<br />
5) 再度、リカバリモードで起動<br />
<br />
<br />
6) 下記のコードを打ち込む(セキュリティを元に戻す)<br />
<br />
<span style="color: magenta;">csrutil enable</span><br />
<br />
<br />
7) 今一度、通常の再起動をする<br />
<br />
<br />
8) Xcodeのインストール、ターミナルに以下を打ち込む<br />
<br />
<span style="color: magenta;">xcode-select --install</span><br />
<br />
<br />
9) homebrewをインストールする<br />
<br />
<span style="color: magenta;">ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"</span><br />
<br />
(せいぜい数分でインストールは終了するはず)<br />
<div>
<br />
<br /></div>
10) XQuartzのインストール。もはやアプリとして付属しなくなった。コードでインストール可能であるようなのでインストール。<br />
<br />
<span style="color: magenta;">brew cask install xquartz</span><br />
<br />
<br />
11) Wineのインストール、下記のコードをターミナルに打ち込む。<br />
<br />
<span style="color: magenta;">brew install wine</span><br />
<br />
<br />
12) 回線速度にもよるが、良好な環境では1時間程度でインストールは終わるだろう。これが終わったら、次のコマンドを打つ。<br />
<br />
<span style="color: magenta;">winecfg</span><br />
<br />
<br />
13) XQuartzのウインドウが表示され、"Gecko package"が必要だから入れてもいいか?と聞かれた。表示されているInstallボタンで承諾するとインストールされる。<br />
<br />
しばらく処理がされた後、XQuartzからWineの環境設定のようなウインドウが表示される。単に一番下のOKをクリックすればよい。<br />
<br />
<br />
14) WinBUGSのインストール<br />
パッチなどを当ててある展開済みのWinBUGSフォルダを用意する。Windows7、Windows8へのインストールも現在はこのやり方で行くしかないことを考えれば、Macでも同様にすればいいだろう。<a href="http://www.wcsmalaysia.org/analysis/WinBUGS_install.htm">適用済みのWinBUGSも公開されている。</a><br />
<br />
下記のコードで不可視フォルダにあるProgram Filesフォルダを開く<br />
<br />
<span style="color: magenta;">open ~/.wine/drive_c/Program\ Files/</span><br />
<br />
ここへWinBUGSフォルダを入れればインストール完了<br />
<div>
<br /></div>
15)RからWinBUGSを実行する。下記のような単純なサンプルコードで試してみる。Wineを経由するので、bugs()内にそのためのコードがたくさん必要。<br />
<br />
<br />
# R2WinBUGSのインストールをお忘れなく<br />
<span style="color: lime;">require(R2WinBUGS)</span><br />
# 真の値は、a=3, b=2, sd=1<br />
<span style="color: lime;">X <- c(1:100)</span><br />
<span style="color: lime;">Y <- rnorm(100, mean=(3 + 2*X), sd=1)</span><br />
<span style="color: lime;">data <- list(X=X, Y=Y)</span><br />
<span style="color: lime;">inits <- function() list(a=0, b=0, tau=1) </span><br />
<span style="color: lime;">parameters <- c("a", "b", "sigma")</span><br />
<span style="color: lime;"><br /></span>
<span style="color: lime;">model <- function() {</span><br />
<span style="color: lime;">a ~ dnorm(0, 1.0E-6)</span><br />
<span style="color: lime;">b ~ dnorm(0, 1.0E-6)</span><br />
<span style="color: lime;">tau ~ dgamma(1.0E-2, 1.0E-2)</span><br />
<span style="color: lime;">for (i in 1:100) {</span><br />
<span style="color: lime;">Y[i] ~ dnorm(mean[i], tau)</span><br />
<span style="color: lime;">mean[i] <- a + b*X[i] }</span><br />
<span style="color: lime;">sigma <- 1/sqrt(tau)</span><br />
<span style="color: lime;">}</span><br />
<span style="color: lime;">modelpath <- file.path(tempdir(), "model.bug")</span><br />
<span style="color: lime;">write.model(model, modelpath)</span><br />
<span style="color: lime;"><br /></span>
<span style="color: lime;">mcmc <- bugs(</span><br />
<span style="color: lime;">data=data, inits=inits, parameters=parameters, model.file=modelpath, </span><br />
<span style="color: lime;">n.chains=3, n.iter=5000, debug=T,</span><br />
<span style="color: lime;">working.directory=NULL, clearWD=T, useWINE=T, newWINE=T,</span><br />
<span style="color: lime;">WINE="/usr/local/bin/wine", WINEPATH="/usr/local/bin/winepath")</span><br />
<span style="color: lime;"><br /></span>
<span style="color: lime;">print(mcmc)</span> # ちゃんと真の値(a=3, b=2, sigma=1)が推定できたかチェックしよう<br />
<br />
# 今回、opt/localではなくusr/localにパスを通すよう変更する必要が出た。以前のWINE、WINEPATHは/opt/local/bin/になっていたが、ここは/usr/local/bin/に変更していることに注意。<br />
<br />
<br />
17)まだR上で下記のエラーコードが出るが、これはこちら(http://ggorjan.blogspot.jp/2008/10/runnning-r2winbugs-on-mac.html)によると害のないエラーコードらしい。要は推定計算さえ無事に行われていればよいだろう。<br />
err:ole:CoGetClassObject class {0003000a-0000-0000-c000-000000000046} not registered<br />
err:ole:CoGetClassObject class {0003000a-0000-0000-c000-000000000046} not registered<br />
err:ole:CoGetClassObject no class object {0003000a-0000-0000-c000-000000000046} could be created for context 0x3<br />
err:ole:CoReleaseMarshalData IMarshal::ReleaseMarshalData failed with error 0x8001011d<br />
<br />
<br />Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-23406985488681570772016-11-20T20:27:00.002+09:002016-12-16T20:38:51.567+09:00潜水から無事に帰還する装備について考える海中での潜水調査・作業をよくやる人の中には、ヒヤッとした経験のある人も多いだろう。私自身もかつてBC(浮力調節)ジャケットの弁のトラブルで強制浮上になってしまったことがあり、それ以降はとくに気になるようになってきた。<br />
<br />
先週、知っている人の中で<a href="https://www.oist.jp/ja/news-center/news/2016/11/18/27851">重大な潜水事故が起こってしまった</a>。まったくもって今更ではあるけれど、ここでは海中または浮上後の海上でトラブルが起こった際に無事に帰還するための装備について考えてみたい。ご意見も歓迎します。<br />
<br />
まず、流されてしまったか何かで、岸や船から遠い地点に浮上してしまった場合、BCに付属している笛を吹く程度では、よほど凪でもなければ気づいてもらえないだろう。他に入手しやすいものとしては、海面上に長く伸びるフロート(レギュレータからのエアを入れて膨らませられる)があるだろう。私も念のため持っているけれど、しかし実際これで十分に目立つのだろうか。<br />
<br />
他の手段…例えばスマホを完全防水で携帯して、浮上したら電話するというのはどうだろうか。場合によってはアカウントを利用して、PCからスマホを探すことも可能な気がするのだけれど。<br />
<br />
スマホに近いものとして、<a href="http://www.855756.com/child/faq02.html">ココセコム</a>というのがあるそうだ。子供やお年寄りに持たせておいて携帯電話網でマップ上で追跡が可能、料金的にもスマホ利用より安そうです。<br />
<br />
無線を使うという手段は電波法上かなり難しいようだが、合法的にダイビングで使用できる製品もあるようだ(<a href="http://www.safedive.jp/?p=343">トランスポンダ SEAKER_1</a>)<br />
船からのダイビングという形態を前提としているようだけれど、通常使用(自分の母船に連絡)と、緊急使用(母船とはぐれている場合に無線を受けられるすべての船に緊急発信が可能)の2通りの発信が可能というのが優れている。距離的には20 kmまで届くようなので頼もしい。<br />
<br />
<b>(2016.12.16 追記)意識があって海上で操作することが前提ですが、その条件ではかなり有効な製品が出たようです(<a href="http://plbjapan.wixsite.com/plb1">個人用救難信号発信機 PLB</a>)。精度100mで全世界対応。機器そのものは今年3月に日本で認可、防水ケースが先月末に発売された模様。無線局として申請・許可を得る必要があるなど取り扱いは注意ではあります。</b><br />
<br />
では、意識がないけれど海面には浮かんでいる場合どうするか。上記の携帯を利用した方法が可能であったら使えそうではあるけれど、岸から遠くへ流されてしまっていたとしたら(携帯電話エリア外)アウトだろう。<br />
登山用にはもう少し手軽な製品が出ているようだ(<a href="http://www.authjapan.com/about-1">ヒトココ</a>)。<br />
こちらは最大距離 1 kmと短めだけれど、親機を安全なところに確保しておいて、わずか20gの子機を携帯していくという形で使用できる模様。完全生活防水とのことだけれど、ダイビングで使うなら何かしらの防水ケースが別途必要だろう。距離の制限はあるけれど、水面に浮いてさえしてくれればいいところは汎用性は高そうだ(意識の有無にかかわらず…)。<br />
<br />
その他、ビーコンという案もいただいた。<br />
<br />
またバイオロギングで海獣にGPS発信機を付けるという話もあるけれど、あれはダイビングの安全管理には使えないのだろうか?<br />
<br />
しかし、もし沈んだままで浮いていなかったら、これらの手段はどうにもならないですね。。海底で作業している間に意識を失ったとしたら、中性浮力でなくエアは抜いた状態のはず。Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-87352960263167293722016-09-02T18:30:00.002+09:002020-08-29T18:07:22.107+09:00(cowplotパッケージ)研究用にスッキリ簡潔にggplotを描画 & 複数パネル化最近、Rでのグラフ作成の標準になりつつある気がするggplot2パッケージですが、デフォルトのテーマは研究用としては装飾過剰なので、自分用にアレンジしたテーマを使っている人も多いと思います。でもそのテーマを図示の度に毎回引っ張ってくるのはとても面倒。<br />
<br />
それから、複数の異なる種類のグラフを組み合わせて描画するときに、gridExtraパッケージを使うというのがありますが、図示するまではいいけれど、保存する時に行・列数の指定ができないなど、こちらもいろいろ厄介でした。<br />
<br />
この両方を一気に解決してくれたのが<a href="https://cran.r-project.org/web/packages/cowplot/vignettes/">cowplotパッケージ</a>、研究目的のスッキリしたグラフを作成するのに特化したようなパッケージです。ggplot2を基本にして拡張したようなものなのですが、使い方はかなり簡単です。以下、irisデータを用いて図示してみます。<br />
<br />
<span style="color: lime;">data(iris) </span># iris呼び出し<br />
<br />
# いずれも要インストール<br />
<span style="color: lime;">require(ggplot2)</span><br />
<span style="color: lime;">require(cowplot)</span><br />
<br />
初めにggplot2とcowplotの両方を呼び出しておけば、あとは普段通りにggplotで図示するだけで論文ライクなスッキリしたグラフが描かれます。<div>→(2020.08.29追記)cowplotの仕様が変更になっており、下記のようなcowplotスタイルのグラフにするには、<span style="color: #04ff00;">theme_set(theme_cowplot())</span>、を通しておく必要があります。または<span style="color: #04ff00;">ggplot()... + theme_cowplot()</span>のように追加してもOKです。<br />
<div>
<br /></div>
<div>
# x=Sepal.Length, y=Sepal.Widthの散布図をSpecies毎に色分けして描く</div>
<div>
<div>
<span style="color: lime;">g1 <- ggplot(iris, aes(x=Sepal.Length, y=Sepal.Width)) + geom_point(aes(fill=Species), pch=21, size=2) </span></div>
<div>
<span style="color: lime;">g1</span> # 図示</div>
<div>
<br /></div>
<div>
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhfb3vxpK6BBykZze6QRQ-AC8b0J0clb7rISJqX9bpH3mRu1ud30D0Z3sNCflWx00P5m3vpXS2WNHlRv4aNVWhyphenhyphen7ffvVeT2LUUjax7lXqq9WSE8k15ayl4lwicu_sWOaFaQpxhyphenhyphenjlMGYWE/s1600/cowsample1.jpg"><img border="0" height="264" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhfb3vxpK6BBykZze6QRQ-AC8b0J0clb7rISJqX9bpH3mRu1ud30D0Z3sNCflWx00P5m3vpXS2WNHlRv4aNVWhyphenhyphen7ffvVeT2LUUjax7lXqq9WSE8k15ayl4lwicu_sWOaFaQpxhyphenhyphenjlMGYWE/s320/cowsample1.jpg" width="320" /></a></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
# 保存も普通にggplot2と同様にやるだけ</div>
<div>
<span style="color: lime;">ggsave("cowsample1.pdf", g1, height=10, width=12, unit="cm")</span></div>
</div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
次は、複数種類のグラフを並べる場合</div>
<div>
<br /></div>
<div>
# 先ほどのg1と新たにg2を用意し、横並びにします。</div>
<div>
<div>
<span style="color: lime;">g2 <- ggplot(iris, aes(x=Petal.Length, y=Petal.Width)) + geom_point(aes(fill=Species), pch=21, size=2)</span></div>
<div>
<span style="color: lime;"><br /></span></div>
# 複数グラフを合成するにはplot_grid関数を使用します<br />
<div>
<span style="color: lime;">g12 <- plot_grid(g1, g2, labels="auto", align="h") </span></div>
<div>
<span style="color: lime;"><br /></span></div>
<div>
# labels="auto":小文字で各パネルにラベリング("AUTO"だと大文字)</div>
<div>
# labels=c("a","b") のように直接入力することも可能</div>
<div>
# align="h":グラフの上下が揃うようにサイズ調整してくれる</div>
<div>
<br /></div>
<div>
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj_kExBLgSEZFaNRx0WEgJYON96u6iseXNQyVzE-E8Ot_muDEl6BUE2sS_IHeJl2IsUqCh7DLiWcCUIEbaaonSwX8SaB_Hr-rwN8xtQgIaaIa6-jX1AEL0qbfWz2JYtnn4c8KGS0x3QRwU/s1600/cowsample12.jpg"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj_kExBLgSEZFaNRx0WEgJYON96u6iseXNQyVzE-E8Ot_muDEl6BUE2sS_IHeJl2IsUqCh7DLiWcCUIEbaaonSwX8SaB_Hr-rwN8xtQgIaaIa6-jX1AEL0qbfWz2JYtnn4c8KGS0x3QRwU/s1600/cowsample12.jpg" /></a></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
# このグラフのように凡例が共通の場合、一方は消してしまいましょうか。</div>
<div>
<span style="color: lime;">g1n <- g1 + theme(legend.position="none")</span> # ggplotのtheme関数と合わせ技が可能、1パネル目の凡例を消去</div>
<div>
<span style="color: lime;">g1n2 <- plot_grid(g1n, g2, labels="auto", rel_widths=c(1, 1.4)</span><span style="color: lime;">, scale=0.9, vjust=0.2</span><span style="color: lime;">) </span></div>
<div>
# rel_widthsで横幅を1:1.4に変更</div>
<div>
# それからデフォルトのままだと、パネル間がくっつきすぎてしまうので、scaleパラメータで調整します(1→0.9)。vjustはラベル位置(a, b)の調整</div>
<div>
<br /></div>
<div>
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhe3W6AQqZuOhCivCXOyk9-5QI5lxvSiGIAsofQor1me8RSx3T45FpltcSkwobiDuM7GuSdYkRF3irpc5BpgEpIAZdyGK0wlnKo_d2AxHuMFHnLuXfB0ZW4TUT-13xVoSWqdH6Hm0vIwtU/s1600/cowsample1n2_1.jpg"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhe3W6AQqZuOhCivCXOyk9-5QI5lxvSiGIAsofQor1me8RSx3T45FpltcSkwobiDuM7GuSdYkRF3irpc5BpgEpIAZdyGK0wlnKo_d2AxHuMFHnLuXfB0ZW4TUT-13xVoSWqdH6Hm0vIwtU/s1600/cowsample1n2_1.jpg" /></a></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
# もう一つ増やして、二段重ねにしてみます(nrow=2で2行になる:gridExtraパッケージのarrageGrob関数と違って、行・列数の指定が可能。grid.arrange関数があるじゃないかと思うかもですが、ファイル保存できないので実用不可)</div>
<div>
<span style="color: lime;">g3 <- ggplot(iris, aes(x=Species, y=Petal.Width)) + geom_boxplot(aes(fill=Species)) + theme(axis.text.x = element_text(angle=90, vjust=0.5))</span><br />
# 追加のg3はx軸のラベルが重なってしまうので縦にする↑(themeを使用)</div>
<div>
<span style="color: lime;">g123 <- plot_grid(g1, g2, g3, labels="auto", nrow=2, scale=0.9, vjust=0.2</span><span style="color: lime;">) </span></div>
<div>
<br /></div>
<div>
</div>
<div>
<div>
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjcLHEcx0Zz9VjP1WlEQYOpIn0QwDzcTIpnsd8bVMj-AK3SRVHr-mudVKqYMZeh2u_f0Wp7DghoN-gjT0Sg7WhsgpvxjXfVhc-WI4_Q4ArAhUMEMGWRgOxAFh_g8XWSFXMEQprZsxsXbtg/s1600/cowsample123_2.jpg"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjcLHEcx0Zz9VjP1WlEQYOpIn0QwDzcTIpnsd8bVMj-AK3SRVHr-mudVKqYMZeh2u_f0Wp7DghoN-gjT0Sg7WhsgpvxjXfVhc-WI4_Q4ArAhUMEMGWRgOxAFh_g8XWSFXMEQprZsxsXbtg/s1600/cowsample123_2.jpg" /></a></div>
<div>
# ちゃんとこの状態で保存も可能です</div>
<div>
<span style="color: lime;">ggsave("cowsample123.pdf", g124, height=15, width=19, unit="cm")</span></div>
</div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<div>
さらに、ggplotらしく要因で分けたグラフと普通のグラフを組み合わせてみます。縦長と普通サイズの組み合わせをしようとしているのですが、この場合、plot_gridを重ねることで実現できます。ラベルをautoのままにすると重なってしまうので注意が必要です。</div>
</div>
<div>
<br /></div>
<div>
# g1, g3で縦に連結(横に持ってくるので、ラベルがb, cになるようにする)<br />
# 幅も揃えておきます(align="v")<br />
<span style="color: lime;">g13 <- plot_grid(g1, g3, labels=c("b","c"), ncol=1, vjust=0.5, scale=0.9, align="v") </span></div>
<div>
<div>
<br />
<div>
# 次にfacet_wrapで種毎のグラフ(ncol=1で1列になる)</div>
<div>
<div>
</div>
</div>
<span style="color: lime;">g4n <- ggplot(iris, aes(x=Sepal.Length, y=Petal.Length)) + geom_point(aes(fill=Species), pch=21, size=2) + facet_wrap(~ Species, ncol=1) + theme(legend.position="none")</span><br />
<br />
# 両者をさらにplot_gridの中に入れます(これによりplot_gridが入れ子になる)。<br />
# 新たな左のグラフのみにラベルをつけるので、c("a","")のように指定しました。<br />
<span style="color: lime;">g134 <- plot_grid(g4n, g13, labels=c("a",""), nrow=1, rel_widths=c(1, 1.4), scale=0.9, vjust=0.1)</span></div>
<div>
<br /></div>
<div>
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiAKxzcSmWSsz4wuhNTOyeME6-ax9wZiLxqNfcc5vWd5NOIQvIUdpTQXA579CMjuNYfLWLHRDiF-PX99AsTBHfDCOA2kUyrZ1VBviXApBuUdNvsSfKNxq4ja3MkwXyrGVHztHP79hawW3A/s1600/cowsample134_6.jpg"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiAKxzcSmWSsz4wuhNTOyeME6-ax9wZiLxqNfcc5vWd5NOIQvIUdpTQXA579CMjuNYfLWLHRDiF-PX99AsTBHfDCOA2kUyrZ1VBviXApBuUdNvsSfKNxq4ja3MkwXyrGVHztHP79hawW3A/s1600/cowsample134_6.jpg" /></a></div>
</div>
</div>
<div>
</div>
<div>
</div>
<div>
<br />
<br /></div>
<div>
<br /></div>
<div>
単体のggplot2と比べて、とても便利。日常的に使いたくなるパッケージです。<br />
<br />
2017.05.19追記:<br />
theme関数で外枠を追加しようとすると上手く行かず。代わりにcowplotパッケージのpanel_border関数を使えば良い模様。<br />
<span style="color: lime;">+ panel_border(colour=1, size=1)</span> # defaultの線色がなぜかグレーなので黒に変更、sizeで線太さを調節</div>
<div>
<br /></div>
</div>Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-55228111214749988422016-08-31T19:40:00.002+09:002020-04-02T11:03:28.934+09:00RのforループをC++で高速化する(Rcppパッケージ)"時間が掛かるからRでforループを使うな"、"applyファミリーを使え(ほかにsapply, lapplyなど)"とよく言われる。かといって無理に使うと難解なコードになる場合もあるし、せっかく実現しても計測してみたら、むしろforループの方が早かったなんてこともあった(今回の動機)。並列化(複数のCPUコアを用いる)するという手も考えたが、大きなー繋がりのタスクを1回やるだけならよいが、入れ子になっていると呼び出す時間の方が律速になったりもする。<br />
<br />
applyファミリーやforeachなどで置換えにくい計算例として、値を毎回更新するようなプログラムがあるだろう。例えば以下のような計算例("0.9* "さえ無ければ<span style="color: lime;">apply(x, 2, cumsum)</span>で一発なのだが、sapplyとapplyを組み合わせて書いたらforより遅くなった…)。<br />
<br />
<span style="color: lime;">NN <- 10^6</span><br />
<span style="color: lime;">x <- matrix(0, NN, 10) </span><br />
# NN回分の結果を格納するための入れ物(1回あたりは長さ10のベクトル)<br />
<span style="color: lime;">system.time( for(r in 2:(nrow(x))) x[r,] <- 0.9*x[r-1,] + rnorm(10) )</span><br />
# ひとつ前の結果に0.9を掛け、正規乱数を足していくマルコフ連鎖<br />
# ユーザ システム 経過 <br />
# 5.457 0.363 5.778<br />
<br />
手元の環境では6秒弱掛かった。<br />
<br />
いろいろ調べて行き着いたのが <span style="color: lime;">Rcpp</span> というパッケージ(要インストール)。RからC++("しーぷらぷら"と呼ぶみたい、なのでcppなのだろう)で作成した関数を呼び出せるというもの。Cで作られているというパッケージは最近よく見かけるようになったが、自分でカスタマイズできるのはすごい。もっとも、C++を書けるならばですが。それでも部分的にだけでも記述できればだいぶ速くすることができるはず。目的のforループ作成を達成するのに丸一日掛かりましたが、最終目的の計算に掛かる予想時間が数ヶ月だったことを考えれば大幅な時間短縮になりました。<br />
<br />
手順として、まず開発環境のセットアップが必要です。Macの場合はXcode(AppStoreからフリーでダウンロード)とX11(ユーティリティ内にあるX11のアイコンをダブルクリックして立ち上げてライセンスにOKすれば大丈夫なはず)。Windows版はすみませんよく分かりません…コマンドプロンプトでC++をコンパイルする必要があると思います。<br />
次に肝心の、C++で記述したコードが必要です。ただし随所に純粋のC++と異なる書き方を要する箇所があります(私はここでつまづきました)。C++で記述したファイルをRの作業ディレクトリに置いて読み込むと、Rの関数として使用できるようになります。<br />
<br />
<br />
<br />
<br />
<br />
まずは最小限のコード例:(以下、C++のコードは紫にしておきます)<br />
<span style="color: lime;">require(Rcpp)</span> # Rcppパッケージを呼び出す(この行はRに直打ち)<br />
<br />
<span style="color: magenta;">#include <Rcpp.h></span><br />
<span style="color: magenta;">using namespace Rcpp;</span><br />
<span style="color: magenta;">// [[Rcpp::export]]</span><br />
<div>
<span style="color: magenta;"><br /></span></div>
<span style="color: magenta;">double testF(double x) {</span><br />
<span style="color: magenta;"> return(x*x);</span><br />
<span style="color: magenta;">}</span><br />
<div>
<br /></div>
紫で書いた6行分をコードエディタなどで書いて test.cpp というファイル名で作業ディレクトリに保存します( .cpp ってC++ファイルの拡張子なのだと思います)。<br />
C++用コードの初めの3行は決まり文句のようなものです。3行目は用途によって変える場合があるようですが、とりあえず忘れてよさそう。<br />
まずRとの大きな違いはすべての変数型を逐一定義すること。この性質は嫌いじゃないです。むしろ勝手に変数型が変わるのはRが嫌われる点でもありますね。<br />
<br />
double は実数の意味、整数なら int 、のように、使う変数の前に半角スペースを置いて書いて定義する必要があります。<br />
<br />
testF は、いま定義しようとしている自作関数の名前です。括弧の中で x という変数が与えられた時に { } の処理を行う変数を作成することを宣言しています。<br />
<br />
return(x*x); これは x*x の結果を返すことを示しています。もうひとつRと大きな違いとして各行の最後に ;(セミコロン)が入ることです(中括弧 { の後ろ以外)。<br />
<br />
<br />
次にR上で以下のようにして、この自作関数を呼び出します(初回の呼び出しには数秒掛かるかもしれない)。<br />
<span style="color: lime;">sourceCpp("test.cpp")</span><br />
<br />
すると、testF 関数が使用できるようになります。例えば…<br />
<span style="color: lime;">testF(1.41421356) </span># 一夜一夜に人見頃<br />
結果はちゃんと 2 が返ってくるはず。<br />
<br />
<br />
ちなみにRの上でC++コードを走らせることも可能です。""で挟んで、code=code.testで指定してやるだけです。<br />
<span style="color: lime;">sourceCpp(code="</span><br />
<span style="color: lime;">#include <Rcpp.h></span><br />
<span style="color: lime;">using namespace Rcpp;</span><br />
<span style="color: lime;">// [[Rcpp::export]]</span><br />
<span style="color: lime;"><br /></span>
<span style="color: lime;">double testF(double x) {</span><br />
<span style="color: lime;"> return(x*x);</span><br />
<span style="color: lime;">} ")</span><br />
<br />
ただし、Rcpp.h以外のファイルを読み込む必要がある場合のやり方が今ひとつ分からず。inlineパッケージのcxxfunctionを使ってできそうなのですが。<br />
<br />
<br />
<br />
<br />
<br />
<br />
次の例では段階を上げて行きます。<br />
<span style="color: lime;">rnorm</span> 関数をC++で書いてみましょう。なお、乱数の発生アルゴリズムとしてRの場合はMersenne Twister法が使用されていますが、C++でやろうとすると一工夫必要でした。これについては、<a href="http://www.sat.t.u-tokyo.ac.jp/~omi/random_variables_generation.html#Prepare_MT">詳しい解説とコードも配布されているサイトがありましたのでリンクを貼っておきます。</a>MT.hファイルと内包されるinit_genrand関数、genrand_real3関数はこちらから取得しました。(追記:他にも、<a href="http://en.cppreference.com/w/cpp/numeric/random">匿名さんから頂いたコメントで乱数発生アルゴリズムのまとめサイトがあるとのことです。</a>)<br />
<br />
r_normal という名前でC++版rnormを作成します。<br />
<span style="color: #6fa8dc;">下準備としてMt.hと名付けたファイルをRの作業ディレクトリに置きます</span>。下記がRcppから読み込む用のC++コードになります(もはや純粋なC++コードではないので、こう呼びます)。<br />
<br />
<span style="color: magenta;">#include <Rcpp.h></span><br />
<span style="color: magenta;">#include "MT.h"</span><br />
<br />
<span style="color: magenta;">using namespace Rcpp;</span><br />
<span style="color: magenta;">// [[Rcpp::export]]</span><br />
<br />
<span style="color: magenta;">NumericVector r_normal( int NN, double mu, double sigma ) { </span><span style="color: magenta;"><span class="Apple-tab-span" style="white-space: pre;"> </span></span><br />
<span style="color: magenta;"> NumericVector x(NN);</span><br />
<span style="color: magenta;"><span class="Apple-tab-span" style="white-space: pre;"> </span>init_genrand((unsigned)time(NULL));</span><br />
<span class="Apple-tab-span" style="white-space: pre;"><span style="color: magenta;"> </span></span><br />
<span style="color: magenta;"><span class="Apple-tab-span" style="white-space: pre;"> </span>for(int i=0; i<NN; i++) {</span><br />
<span style="color: magenta;"><span class="Apple-tab-span" style="white-space: pre;"> </span>x[i] = mu + sigma*sqrt( -2.0*log(genrand_real3()) ) * sin( 2.0*M_PI*genrand_real3() );</span><br />
<span style="color: magenta;"><span class="Apple-tab-span" style="white-space: pre;"> </span>}</span><br />
<span style="color: magenta;"><span class="Apple-tab-span" style="white-space: pre;"> </span>return(x);</span><br />
<span style="color: magenta;">}</span><br />
<div>
<br /></div>
新たに加わった include "MT.h" で上記のファイルを詠みこみます<br />
<br />
<span style="color: magenta;">#include <stdio.h></span><br />
<span style="color: magenta;">#include "MT.h"</span><br />
<br />
次に、この行ですが、これがクセモノです。<br />
<span style="color: magenta;">NumericVector r_normal( int NN, double mu, double sigma ) {</span><br />
<br />
NumericVector は返り値がベクトルになるように定義しますという意味です(Rcpp特有の表現:cf. <a href="https://cran.r-project.org/web/packages/Rcpp/vignettes/Rcpp-quickref.pdf">Rcpp Quick reference</a>)。そして r_normal という名前の自作関数を作ろうとします。<br />
<br />
int NN, double mu, double sigma は、ベクトルの要素数 NN(整数)、平均が mu(実数)、標準偏差が sigma(実数)の値を用いることを宣言しています。<br />
<br />
<span style="color: magenta;">NumericVector x(NN);</span><br />
次にこれですが、xという、要素数が NN のベクトル(初期値は全部0になっている)を用意することを意味しています(Rcpp特有の定義のようだ)。<br />
<br />
<span style="color: magenta;">init_genrand((unsigned)time(NULL));</span><br />
これは乱数のタネの指定、現在時刻を秒精度で使用します。<br />
<br />
<span style="color: magenta;"><span class="Apple-tab-span" style="white-space: pre;"> </span>for(int i=0; i<NN; i++) {</span><br />
<div>
C++版のforループ。整数 i を 0からNN未満まで適用しますという意味。0から始めるのが基本である模様(誤解があるかもしれないけれど、とりあえずこのコードではOK)。</div>
<div>
<br /></div>
<div>
<span style="color: #674ea7;"><span class="Apple-tab-span" style="color: #674ea7; white-space: pre;"> </span><span style="color: magenta;">x[i] = mu + sigma*sqrt( -2.0*log(genrand_real3()) ) * sin( 2.0*M_PI*genrand_real3() );</span></span></div>
<div>
<span style="color: #674ea7;"><br /></span></div>
xのi番目の要素に正規乱数を格納するコード。右辺はごちゃごちゃしていますが、genrand_real3が一様乱数を返す関数です。詳しい説明が<a href="http://www.sat.t.u-tokyo.ac.jp/~omi/random_variables_generation.html#Gauss">参照元</a>にあります。<br />
return(x); で x を返り値で返すようにすれば、r_normal関数が使えるようになります。<br />
<br />
rnormと速度比較をしてみます。<br />
<br />
<span style="color: lime;">NN <- 10^7</span><br />
<span style="color: lime;">system.time(rnorm(n, 0, 1))</span><br />
# ユーザ システム 経過 <br />
# 0.677 0.006 0.679<br />
<div>
<br /></div>
<span style="color: lime;">system.time(r_normal(n, 0, 1))</span><br />
# ユーザ システム 経過 <br />
# 0.363 0.005 0.367<br />
<br />
速度は2倍弱といったところ。もう少し期待したかったですが、それでもベクトル化されているRコードよりもさらに高速です。<br />
<br />
<br />
一方、Rの関数をC++側に渡すことも可能です。Rの関数をC++で計算させることができます。やり方としては、cppFunction関数を用いてrnormのC++版、rnormCppを定義します。<br />
<span style="color: lime;">cppFunction( "</span><br />
<span style="color: lime;">NumericVector rnormCpp(int N, double me, double sd) { </span><br />
<span style="color: lime;">return(rnorm(N, me, sd)); </span><br />
<span style="color: lime;">}" )</span><br />
<br />
定義の仕方などはだいたい同じですが、rnormが使われていることに注意です。これはC++の関数ではなく、たしかにRのrnorm関数です。<br />
これを通すと、rnormCppという関数が使えるようになります。気になる速度ですが…<br />
<br />
<span style="color: lime;">system.time(rnormCpp(10^7, 0, 1))</span><br />
# ユーザ システム 経過 <br />
# 0.577 0.019 0.594<br />
<br />
Rのrnormより少し速くはなりましたが、C++で1から作るより遅いです。それでも手軽に計算速度を向上させることができるので、何かと試したくなる技です。<br />
<br />
<div>
<br />
<br />
<br /></div>
<br />
ではいよいよ、本題のforループを r_markov 関数として記述してみます(同様にファイル保存します)。<br />
<span style="color: magenta;"><br /></span>
<span style="color: magenta;">#include <Rcpp.h></span><br />
<span style="color: magenta;">#include "MT.h"</span><br />
<span style="color: magenta;"><br /> using namespace Rcpp;<br /> // [[Rcpp::export]]<br /><br /> NumericMatrix r_markov( int repl, int length, double mu, double sigma ) {</span><br />
<span style="color: magenta;"></span><span class="Apple-tab-span" style="color: magenta; white-space: pre;"> </span><span style="color: magenta;">NumericMatrix x(repl, length);</span><br />
<span style="color: magenta;"><span class="Apple-tab-span" style="white-space: pre;"> </span>int i,j;</span><br />
<span class="Apple-tab-span" style="white-space: pre;"><span style="color: magenta;"> </span></span><br />
<span style="color: magenta;"><span class="Apple-tab-span" style="white-space: pre;"> </span>init_genrand((unsigned)time(NULL));</span><br />
<span class="Apple-tab-span" style="white-space: pre;"><span style="color: magenta;"> </span></span><br />
<span style="color: magenta;"><span class="Apple-tab-span" style="white-space: pre;"> </span>for(i=1; i<repl; i++) {</span><br />
<span style="color: magenta;"><span class="Apple-tab-span" style="white-space: pre;"> </span>for(j=0; j<length; j++) {</span><br />
<span style="color: magenta;"><span class="Apple-tab-span" style="white-space: pre;"> </span>x(i,j) = 0.9*x(i-1,j) + mu + sigma*sqrt( -2.0*log(genrand_real3()) ) * sin( 2.0*M_PI*genrand_real3() );</span><br />
<span style="color: magenta;"><span class="Apple-tab-span" style="white-space: pre;"> </span>}</span><br />
<span style="color: magenta;"><span class="Apple-tab-span" style="white-space: pre;"> </span>}</span><br />
<span style="color: magenta;"><span class="Apple-tab-span" style="white-space: pre;"> </span>return(x);</span><br />
<span style="color: magenta;">}</span><br />
<br />
繰り返しを各行に割り当てるため、int repl、各回のベクトル要素数は int length です。<br />
ひとつ前のr_normal関数との違いは、まず毎回の結果を行列に格納している点です。<br />
NumericMatrixが行列型を示しています。<br />
一行目でもxをrepl行、length列の行列として定義します。NumericMatrix x(repl, length);<br />
毎回の繰り返しでは、ひとつ前の結果に0.9を掛けて正規乱数を足していくので、iのforループは1から始まるようにしています(1行目は0のままで放置、Rコードで除去します)。<br />
もう一つ、分かりにくい違いが添字の書き方です。さっきはx[i] = ...と肩のある括弧を使いましたが、行列の場合はx(i,j)と普通の括弧を使用するようです。C++の解説を見るとx[i][j]とありますが、こうするとエラーになりました。x[i,j]もダメです。<br />
<br />
<br />
<br />
それでは冒頭のR版forループと速度比較してみます。<br />
<span style="color: lime;">NN <- 10^6</span><br />
<span style="color: lime;">system.time(r_markov(NN+1, 10, 0, 1)[-1,])</span><br />
# ユーザ システム 経過 <br />
# 0.439 0.024 0.460<br />
<br />
R版では5.778秒だったので10倍以上早いです。C++の速さを思い知らされました。<br />
<br />
Rcpp、少しずつ部分的にC++化できるので、練習にも持って来いだと思いました。たぶん純粋C++に移行することはなくハイブリッドでやっていくことになりそうな気がします。<br />
<br />Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com5tag:blogger.com,1999:blog-6412094354437530509.post-72208014844594708032015-12-11T15:27:00.001+09:002016-09-19T18:39:23.704+09:00MacOSX10.11 (El Capitan) でWinBUGSを動かす (Wine, R2WinBUGS使用)<a href="http://nhkuma.blogspot.jp/2012/12/macosx108-mountain-lion-winbugs-wine.html">OSX10.8以降、長らくアップデートしていなかったので</a>、更新することにしました。もういい加減WinBUGSは卒業したくはあるのですが、依然StanやJAGSでは通らないコードも残っておりまだ手放せずにいます。早くStanが離散パラメータを直接扱えるようになってほしい…<br />
<br />
クリーンインストールのEl Capitanに入れることを想定して手順を挙げておきます。まず、X11とHomebrew(パッケージ化されていないアプリを容易にインストールするための補助ツール?)を入れ、Homebrewを用いてWine(非Windows OS上でWin専用アプリケーションを実行する環境)をインストール、Wineのディレクトリ内にWinBUGSをインストールするという流れです。以前はMacPortsを使用していましたが、Homebrewの方がはるかに簡単のようです。以前利用できていたユーザの上の階層がEl Capitanでは使いにくくなったというのも大きな理由です。<br />
<br />
なお、インストール作業はターミナルからUNIXコマンドを打ちながらのもの。sudoなどのコマンドは注意深く扱う必要があるようなので、チャレンジする際には慎重に。またXcodeもソフトウェア開発に使うような類のツールなので取り扱い注意です。参考にする際は、この辺りを理解の上、自己責任でお願いします…。<br />
<br />
(下記、Rコードは緑、ターミナルのコードは紫にしてみます)<br />
<div>
<br />
cf. Windows10でWinBUGSを使用するには一手間必要。Program Filesが変更不能になったため、User以下のフォルダにWinBUGSを入れるしか無い。<br />
<span style="color: lime;">bugs(..., bugs.directory="C:/Users/ゆーざぁ名/Documents/WinBUGS14/")</span><br />
のようにして、WinBUGS14.exeファイルを置いている階層を指定してやる必要がある(例は、ドキュメントフォルダ内にWinBUGS14フォルダを置いて、その直下のWinBUGS14.exeを呼び出す場合。"ゆーざぁ名"は御自身の使用しているものに置き換えてください)。</div>
<br />
cf. usr/local階層に変更を加えるための認証は、OSのマイナーアップグレードの度に必要になるかもしれない。その場合、ターミナルで以下のコードを打ち込む<br />
<span style="color: magenta;">sudo chown $(whoami):admin /usr/local && sudo chown -R $(whoami):admin /usr/local</span><br />
<br />
<br />
<br />
**********************<br />
0)実行環境<br />
・OSX10.11 (El Capitan) 搭載のMac<br />
<br />
1)App StoreのApple IDを設定しておく<br />
<br />
<br />
2)「アプリケーション」→「ユーティリティ」にある「ターミナル」を立ち上げる。<br />
すると、冒頭にこのように出ている。<br />
<br />
コンピュータ名:~ ユーザ名$<br />
<br />
このドルマーク $ の後にコマンドを打っていく。<br />
なお、インストールに関わるところでパスワードを求められるが、その都度、自分のアカウントのパスワードを入れる。<br />
(以降、パスワードを入れる作業は説明を省略)<br />
<br />
<br />
3) 下準備、/usr/local/フォルダを作る、ロックをいったん外して操作をするという動作のようなので、推奨されていない動作だということをお忘れなく。<br />
<br />
リカバリモード(⌘+R を押しながら起動)で起動し、ターミナルを立ち上げる<br />
<br />
<br />
4) ターミナルの $ マークの直後に、下記のコードを打ち込む(コピペでOK、以降も同様)<br />
<div>
<br /></div>
<div>
<span style="color: magenta;">csrutil disable</span></div>
<div>
<br /></div>
<br />
5) 通常の再起動をする<br />
<br />
<br />
6) ターミナルに下記を打ち込む(改行されて見えているだろうが、改行無しで打ち込む)<br />
<br />
<span style="color: magenta;">sudo mkdir /usr/local && sudo chflags norestricted /usr/local && sudo chown $(whoami):admin /usr/local && sudo chown -R $(whoami):admin /usr/local</span><br />
<br />
<br />
7) 再度、リカバリモードで起動<br />
<br />
<br />
8) 下記のコードを打ち込む<br />
<br />
<span style="color: magenta;">csrutil enable</span><br />
<br />
<br />
9) 今一度、通常の再起動をする<br />
<br />
<br />
10) homebrewをインストールする<br />
<br />
<span style="color: magenta;">ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"</span><br />
<br />
(せいぜい数分でインストールは終了するはず)<br />
<br />
<br />
11) Wineのインストールを試みる、下記のコードをターミナルに打ち込む<br />
<br />
<span style="color: magenta;">brew install wine</span><br />
<br />
数分ほどでエラーメッセージが出る、メッセージの最後に下記のように書かれている<br />
To continue, you must install Xcode from the App Store, or the CLT by running:<br />
xcode-select --install<br />
<br />
先にXcodeをインストールしておけば回避できるのだろうが、このやり方でインストールするほうがむしろ手間が省けるだろう。<br />
<br />
<br />
12) 上記のメッセージ通り、下記のコードを打ち込む<br />
<span style="color: magenta;">xcode-select --install</span><br />
<br />
すると、AppStoreからXcodeをインストールしてもよいかと聞かれるのでOKをする。<br />
<br />
<br />
13) インストールが終わったら、再びWineのインストールを試みる<br />
<br />
<span style="color: magenta;">brew install wine</span><br />
<br />
<br />
14) 回線速度にもよるが、良好な環境では30分程度でインストールは終わるだろう。これが終わったら、次のコマンドを打つ。<br />
<br />
<span style="color: magenta;">winecfg</span><br />
<br />
しばらく処理音が聞こえた後、X11からWineの環境設定のようなウインドウが表示される。単に一番下のOKをクリックすればよい。<br />
ターミナルにはエラーメッセージがいくつか出ているが気にしなくてよさそうだ。<br />
<div>
Wineのインストールはたったこれだけで終わり…MacPortsの時の面倒を思えば、Wine自体のインストールはずっと容易になった。</div>
<div>
<br /></div>
<div>
<br /></div>
15) WinBUGSのインストール<br />
パッチなどを当ててある展開済みのWinBUGSフォルダを用意する。Windows7、Windows8へのインストールも現在はこのやり方で行くしかないことを考えれば、Macでも同様にすればいいだろう。<a href="http://www.wcsmalaysia.org/analysis/WinBUGS_install.htm">適用済みのWinBUGSも公開されている。</a><br />
<br />
下記のコードで不可視フォルダにあるProgram Filesフォルダを開く<br />
<br />
<span style="color: magenta;">open ~/.wine/drive_c/Program\ Files/</span><br />
<br />
ここへWinBUGSフォルダを入れればインストール完了<br />
<br />
<br />
16)RからWinBUGSを実行する。下記のような単純なサンプルコードで試してみる。Wineを経由するので、bugs()内にそのためのコードがたくさん必要。<br />
<br />
# R2WinBUGSのインストールをお忘れなく<br />
<span style="color: lime;">require(R2WinBUGS)</span><br />
# 真の値は、a=3, b=2, sd=1<span style="color: lime;"><br /></span>
<span style="color: lime;">X <- c(1:100)</span><br />
<span style="color: lime;">Y <- rnorm(100, mean=(3 + 2*X), sd=1)</span><br />
<span style="color: lime;">data <- list(X=X, Y=Y)</span><br />
<span style="color: lime;">inits <- function() list(a=0, b=0, tau=1) </span><br />
<span style="color: lime;">parameters <- c("a", "b", "sigma")</span><br />
<span style="color: lime;"><br /></span>
<span style="color: lime;">model <- function() {</span><br />
<span style="color: lime;">a ~ dnorm(0, 1.0E-6)</span><br />
<span style="color: lime;">b ~ dnorm(0, 1.0E-6)</span><br />
<span style="color: lime;">tau ~ dgamma(1.0E-2, 1.0E-2)</span><br />
<span style="color: lime;">for (i in 1:100) {</span><br />
<span style="color: lime;">Y[i] ~ dnorm(mean[i], tau)</span><br />
<span style="color: lime;">mean[i] <- a + b*X[i] }</span><br />
<span style="color: lime;">sigma <- 1/sqrt(tau)</span><br />
<span style="color: lime;">}</span><br />
<span style="color: lime;">modelpath <- file.path(tempdir(), "model.bug")</span><br />
<span style="color: lime;">write.model(model, modelpath)</span><br />
<span style="color: lime;"><br /></span>
<span style="color: lime;">mcmc <- bugs(</span><br />
<span style="color: lime;">data=data, inits=inits, parameters=parameters, model.file=modelpath, </span><br />
<span style="color: lime;">n.chains=3, n.iter=5000, debug=T,</span><br />
<span style="color: lime;">working.directory=NULL, clearWD=T, useWINE=T, newWINE=T,</span><br />
<span style="color: lime;">WINE="/usr/local/bin/wine", WINEPATH="/usr/local/bin/winepath")</span><br />
<span style="color: lime;"><br /></span>
<span style="color: lime;">print(mcmc) </span># ちゃんと真の値(a=3, b=2, sigma=1)が推定できたかチェックしよう<br />
<br />
# 今回、opt/localではなくusr/localにパスを通すよう変更する必要が出た。以前のWINE、WINEPATHは/opt/local/bin/になっていたが、ここは/usr/local/bin/に変更していることに注意。<br />
<br />
<br />
17)まだR上で下記のエラーコードが出るが、これはこちら(http://ggorjan.blogspot.jp/2008/10/runnning-r2winbugs-on-mac.html)によると害のないエラーコードらしい。要は推定計算さえ無事に行われていればよいだろう。<br />
err:ole:CoGetClassObject class {0003000a-0000-0000-c000-000000000046} not registered<br />
err:ole:CoGetClassObject class {0003000a-0000-0000-c000-000000000046} not registered<br />
err:ole:CoGetClassObject no class object {0003000a-0000-0000-c000-000000000046} could be created for context 0x3<br />
<div>
err:ole:CoReleaseMarshalData IMarshal::ReleaseMarshalData failed with error 0x8001011d</div>
Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-31543405076518745432015-01-29T17:35:00.000+09:002015-01-29T19:17:24.780+09:00RでGIS その 1:シェープファイル操作、図示RでのGIS操作、いずれまとめようと思いつつ放ったらかしてました。だれにでも有用そうなものから少しずつアップしていく予定です。<br />
<br />
基本の関数の備忘録、とくにshapefileの読み込みと書き出しの関数が長くて忘れてしまいがちです。。<br />
<div>
<br /></div>
<div>
<span style="color: lime;">require(maptools) </span># shapefileの読み込みなどに用いるパッケージ(要インストール)</div>
<div>
<span style="color: lime;">shape <- readShapeSpatial(file.choose())</span> # 読み込み & .shpファイルをメニューで選択</div>
<div>
<br /></div>
<div>
<span style="color: lime;">plot(shape) </span># 図示もできます</div>
<div>
<span style="color: lime;">zoom(shape)</span> # 自分で選んだ範囲を拡大する場合。</div>
<div>
# このコマンドを打った後に、図のウィンドウ上で対角線の端と端をクリックで選択するという、Rらしからぬ操作法でズームします</div>
<div>
<span style="color: lime;">plot(shape2, add=T)</span> # 他のファイルshape2を重ねて図示したい場合</div>
<div>
<br /></div>
<div>
<span style="color: lime;">str(shape, 5) </span># shapeの中身を眺める場合("5"くらいに制限しておかないとコンソールが溢れて大変なことになる)</div>
<div>
<span style="color: lime;">shape@data </span># shapeファイルのデータを取り出す場合(@dataの中身はデータフレーム)</div>
<div>
<br /></div>
<div>
<br /></div>
<div>
# 手持ちのデータフレームDataからshapefileを作る場合(LongitudeとLatitudeの列を含むデータとします)</div>
<div>
# GPSデータはWGS84(133.33333のような表記)にするのが原則です</div>
<div>
<span style="color: lime;">require(sp) </span># maptoolsを使用していれば、新たに呼び出さなくてよいはず</div>
<div>
<span style="color: lime;">coordinates(Data) <- c("Longitude", "Latitude")</span> # このようにGPS列を指定すると空間データ化する</div>
<div>
# 変な感じがするかもしれないですが、x, yの順番なのでLongitudeを先に書きます</div>
<div>
<br /></div>
<div>
<span style="color: lime;">bbox(Data) </span># cf. これをやるとデータの四隅(最少・最大)が分かります</div>
<div>
<br /></div>
<div>
#もしデータがグリッド状に揃っている場合はグリッドに変換することができる</div>
<div>
<span style="color: lime;">gridded(Data) <- TRUE</span> # そうでない場合はエラーになるはず</div>
<div>
<br /></div>
<div>
# データの書き出し</div>
<div>
<span style="color: lime;">writeSpatialShape(Data, "</span>ファイル名<span style="color: lime;">.shp")</span> # readの場合とはShapeとSpatialの順番が逆!</div>
<div>
<br /></div>
<div>
<br /></div>
Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-32491367701085868392014-06-10T18:20:00.003+09:002015-12-04T10:09:39.355+09:00累積ロジットとGLM二項分布の比較・再&続(うっかり、同様の検証記事を消してしまったので、ついでにアップデートします)<br />
<a href="http://nhkuma.blogspot.jp/2014/06/r-ordinal.html">前回に引き続き</a>、段階的なカテゴリーデータのモデリングに用いられる累積ロジット(cumulative logit)、<br />
例えば、悪い、ふつう、よい、のようなデータを関連しそうな要因で解析するような場合に用いる。<br />
<br />
しかし、悪い=0、ふつう=1、よい=2、のように数値化してしまえば、二項分布型のGLMではダメなのだろうか?たぶん、間隔のいびつなデータでは累積ロジットの方が適しているのだろうが。<br />
<br />
テストデータを用いて、両者の推定と推定値の求め方を比較する。<br />
<br />
# まず解析用データの設定<br />
<span style="color: lime;">logistic <- function(xx) 1 / (1 + exp(-xx))</span><br />
<span style="color: lime;">N <- 3 </span># 最大3、つまり 0, 1, 2, 3 の値を取りうる<br />
<span style="color: lime;">X <- rep(c(1:10), each=10)</span><br />
<span style="color: lime;">Y1 <- rbinom(100, N, prob=logistic(-5 + 0.8*X)) </span># logisticの中身を基にして、二項乱数を発生<br />
<span style="color: lime;">Y2 <- factor(Y1, ordered=T) </span># 累積ロジット用に、ランク化した応答変数を用意<br />
<span style="color: lime;">D <- data.frame(X, Y1, Y2)</span><br />
<div>
<br /></div>
<div>
# 解析の実行</div>
<span style="color: lime;">require(ordinal)</span><br />
<span style="color: lime;">require(VGAM)</span><br />
<span style="color: lime;">M1 <- glm(cbind(Y1, N-Y1) ~ X, family=binomial, data=D) </span># 二項分布のGLM<br />
<span style="color: lime;">M2 <- clm(Y2 ~ X, data=D) </span># 累積ロジット<br />
<span style="color: lime;">M3 <- vglm(Y2 ~ X, family=cumulative(parallel=T), data=D)</span><br />
# 比較用にvglm版、parallel=Tで切片のみ複数になる、Fにすると回帰係数もランクごとに推定(切片のみランク毎の場合、比例オッズモデルともいうようだ)<br />
<div>
<br /></div>
<div>
### GLM</div>
<div>
<div>
<span style="color: lime;">summary(M1)</span></div>
<div>
# Coefficients:</div>
<div>
# Estimate Std. Error z value Pr(>|z|) </div>
<div>
# (Intercept) -5.9324 0.6644 -8.929 <2e-16 ***</div>
<div>
# X 0.9241 0.1011 9.141 <2e-16 ***</div>
<div>
# ---</div>
<div>
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1</div>
<div>
#</div>
<div>
# (Dispersion parameter for binomial family taken to be 1)</div>
<div>
#</div>
<div>
# Null deviance: 287.723 on 99 degrees of freedom</div>
<div>
# Residual deviance: 84.861 on 98 degrees of freedom</div>
<div>
# AIC: 139.14</div>
</div>
<div>
<br /></div>
<div>
### 累積ロジット(clm)</div>
<div>
<div>
<span style="color: lime;">summary(M2)</span></div>
<div>
# link threshold nobs logLik AIC niter max.grad cond.H </div>
<div>
# logit flexible 100 -67.06 142.12 5(0) 2.20e-08 2.6e+03</div>
<div>
#</div>
<div>
# Coefficients:</div>
<div>
# Estimate Std. Error z value Pr(>|z|) </div>
<div>
# X 1.2295 0.1659 7.41 1.26e-13 ***</div>
<div>
# ---</div>
<div>
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1</div>
<div>
#</div>
<div>
# Threshold coefficients:</div>
<div>
# Estimate Std. Error z value</div>
<div>
# 0|1 5.955 0.913 6.522</div>
<div>
# 1|2 8.082 1.130 7.153</div>
<div>
# 2|3 9.716 1.283 7.575</div>
</div>
<div>
<br /></div>
<div>
<div>
### 累積ロジット(vglm)<span style="color: #6fa8dc;"># 回帰係数の正負が逆になる</span></div>
<div>
<span style="color: lime;">summary(M3)</span></div>
<div>
# Coefficients:</div>
<div>
# Estimate Std. Error z value</div>
<div>
# (Intercept):1 5.9548 0.92423 6.4430</div>
<div>
# (Intercept):2 8.0821 1.14814 7.0393</div>
<div>
# (Intercept):3 9.7158 1.31086 7.4117</div>
<div>
# X -1.2295 0.16832 -7.3046 </div>
<div>
#</div>
<div>
# Number of linear predictors: 3 </div>
<div>
#</div>
<div>
# Names of linear predictors: logit(P[Y<=1]), logit(P[Y<=2]), logit(P[Y<=3])</div>
<div>
#</div>
<div>
# Dispersion Parameter for cumulative family: 1</div>
<div>
#</div>
<div>
# Residual deviance: 134.1208 on 296 degrees of freedom <span style="color: #3d85c6;"># 過小分散気味、clmには無い情報!</span></div>
<div>
#</div>
<div>
# Log-likelihood: -67.06038 on 296 degrees of freedom</div>
</div>
<div>
<br /></div>
<div>
<br /></div>
<div>
### 推定値を得るには少し手間がかかる</div>
<div>
<div>
<span style="color: lime;">pre1 <- round(3*fitted(M1))</span> # GLM:試行回数*確率、を整数値に丸める</div>
<div>
<span style="color: lime;">pre2 <- as.numeric(predict(M2,type="class")$fit) - 1 </span></div>
<div>
# 累積ロジット(clm):predict関数で推定ランクを求め、</div>
<div>
# ラベルに合うように 1 を引く(ランクは1,2,3,4、元の値は0,1,2,3なので)</div>
<div>
<span style="color: lime;">pre3 <- apply(fitted(M3), 1, which.max) - 1</span></div>
<div>
# 累積ロジット(vglm):何番目のランクの確率が最大かを求め、</div>
<div>
# ラベルに合うように 1 を引く</div>
<div>
<br /></div>
<div>
<span style="color: lime;">cbind(X, Y1, pre1, pre2, pre3)</span> </div>
<div>
# 推定は完全ではないが、GLMと累積ロジットの推定結果は近いものになった。</div>
</div>
<div>
<div>
X Y1 pre1 pre2 pre3</div>
<div>
1 1 0 0 0 0</div>
<div>
2 1 0 0 0 0</div>
<div>
3 1 0 0 0 0</div>
<div>
4 1 0 0 0 0</div>
<div>
5 1 0 0 0 0</div>
<div>
(中略)</div>
<div>
51 6 1 1 1 1</div>
<div>
52 6 1 1 1 1</div>
<div>
53 6 2 1 1 1</div>
<div>
54 6 2 1 1 1</div>
<div>
55 6 2 1 1 1</div>
<div>
56 6 1 1 1 1</div>
<div>
57 6 2 1 1 1</div>
<div>
58 6 1 1 1 1</div>
<div>
59 6 1 1 1 1</div>
<div>
60 6 0 1 1 1</div>
<div>
61 7 3 2 2 2</div>
<div>
62 7 1 2 2 2</div>
<div>
63 7 1 2 2 2</div>
<div>
64 7 3 2 2 2</div>
<div>
65 7 3 2 2 2</div>
<div>
66 7 1 2 2 2</div>
<div>
67 7 3 2 2 2</div>
<div>
68 7 2 2 2 2</div>
<div>
69 7 1 2 2 2</div>
<div>
70 7 3 2 2 2</div>
<div>
71 8 3 2 3 3</div>
<div>
72 8 2 2 3 3</div>
<div>
73 8 3 2 3 3</div>
<div>
74 8 3 2 3 3</div>
<div>
75 8 2 2 3 3</div>
<div>
76 8 1 2 3 3</div>
<div>
77 8 1 2 3 3</div>
<div>
78 8 3 2 3 3</div>
<div>
79 8 3 2 3 3</div>
<div>
80 8 3 2 3 3</div>
<div>
81 9 2 3 3 3</div>
<div>
82 9 1 3 3 3</div>
<div>
83 9 3 3 3 3</div>
<div>
(後略)</div>
</div>
<div>
<br /></div>
<div>
<br /></div>
<div>
### では、<span style="color: #3d85c6;">推定確率の曲線を描くにはどうしたらよいか?</span></div>
<div>
# これがけっこう面倒くさい。まずはvglmを用いてチェックする。</div>
<div>
<br /></div>
<div>
# vglmをfitted()すると、ランク毎の確率が出てくる</div>
<div>
<span style="color: lime;">head(fitted(M3))</span> # head()で頭出しする</div>
<div>
0 1 2 3</div>
<div>
<div>
1 0.991209875 0.007734525 0.0008493615 0.0002062383</div>
<div>
2 0.991209875 0.007734525 0.0008493615 0.0002062383</div>
<div>
3 0.991209875 0.007734525 0.0008493615 0.0002062383</div>
<div>
4 0.991209875 0.007734525 0.0008493615 0.0002062383</div>
<div>
5 0.991209875 0.007734525 0.0008493615 0.0002062383</div>
<div>
6 0.991209875 0.007734525 0.0008493615 0.0002062383</div>
</div>
<div>
<br />
# clmでは、<br />
<span style="color: lime;">M2@fitted.values</span> で同様にランク毎の確率が得られる。<br />
<br />
<br /></div>
<div>
<br /></div>
<div>
# まず、この出来合いの推定値を図示してみる</div>
<div>
<div>
<span style="color: lime;">plot(X, fitted(M3)[,1], xlim=c(0,10), ylim=c(0,1), col="lightblue", ylab="probability") </span># 0</div>
<div>
<span style="color: lime;">points(X, fitted(M3)[,2], col="turquoise")</span> # 1</div>
<div>
<span style="color: lime;">points(X, fitted(M3)[,3], col="royalblue") </span># 2</div>
<div>
<span style="color: lime;">points(X, fitted(M3)[,4], col="darkblue") </span># 3</div>
<div>
<br /></div>
<div>
# モデルの構造を考えると推定確率曲線は次の計算でいいはず</div>
<div>
#(coef(M3)[4]は回帰係数だが、正負が逆なので - を付ける)</div>
<div>
# 累積ロジットの名前通り、累積で表されている</div>
<div>
<span style="color: lime;">curve(1 - logistic(-coef(M3)[4]*x - coef(M3)[1]), </span><br />
<span style="color: lime;"> add=T, lwd=2, col="lightblue")</span> # ランク0: 1 - ランク1の累積確率</div>
<div>
<span style="color: lime;">curve(logistic(-coef(M3)[4]*x - coef(M3)[1]) - logistic(-coef(M3)[4]*x - coef(M3)[2]), </span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>add=T, lwd=2, col="turquoise")</span> # ランク1: ランク1の累積確率 - ランク2の累積確率</div>
<div>
<span style="color: lime;">curve(logistic(-coef(M3)[4]*x - coef(M3)[2]) - logistic(-coef(M3)[4]*x - coef(M3)[3]),</span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span> add=T, lwd=2, col="royalblue")</span> # ランク2: ランク2の累積確率 - ランク3の確率</div>
<div>
<span style="color: lime;">curve(logistic(-coef(M3)[4]*x - coef(M3)[3]), </span><br />
<span style="color: lime;"> add=T, lwd=2, col="darkblue")</span> # ランク3: ランク3の確率(もはや累積でない)<br />
<br />
<span style="color: lime;">curve(logistic(coef(M1)[1] + coef(M1)[2]*x), add=T, col="tomato", lwd=2)</span> # 比較用にGLMも<br />
<br /></div>
</div>
<div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhwSSPOdd9B3NbN9rMMEBvQub_jD7Fym8G9I_WeCd1MRSja10rGXOkfdMBkYhTZdN0Q__rW4azk42QVBCmA8iI-Nv3MAleE8TPY4hXxq32fUR8XM2-RiSc8p3yiHDtztxOk2w-PaIUW8rc/s1600/%E3%82%B9%E3%82%AF%E3%83%AA%E3%83%BC%E3%83%B3%E3%82%B7%E3%83%A7%E3%83%83%E3%83%88+2014-06-10+18.51.27.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="299" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhwSSPOdd9B3NbN9rMMEBvQub_jD7Fym8G9I_WeCd1MRSja10rGXOkfdMBkYhTZdN0Q__rW4azk42QVBCmA8iI-Nv3MAleE8TPY4hXxq32fUR8XM2-RiSc8p3yiHDtztxOk2w-PaIUW8rc/s1600/%E3%82%B9%E3%82%AF%E3%83%AA%E3%83%BC%E3%83%B3%E3%82%B7%E3%83%A7%E3%83%83%E3%83%88+2014-06-10+18.51.27.png" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
# 色が薄い〜濃いにかけて、それぞれ、0, 1, 2, 3 になる確率、赤はGLMの確率</div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
# ちゃんと関数による推定プロットと、自前で作った推定曲線が一致した。計算の仕方はこれでよさそうだ。</div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
</div>
<div>
<div>
# ちなみにclmの場合はこう計算する(回帰係数の前の - が不要)</div>
<div>
<span style="color: lime;">curve(1 - logistic(coef(M2)[4]*x - coef(M2)[1]), add=T, lwd=2, col="lightblue") </span># 0</div>
<div>
<span style="color: lime;">curve(logistic(coef(M2)[4]*x - coef(M2)[1]) - logistic(coef(M2)[4]*x - coef(M2)[2]), </span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>add=T, lwd=2, col="turquoise") </span># 1</div>
<div>
<span style="color: lime;">curve(logistic(coef(M2)[4]*x - coef(M2)[2]) - logistic(coef(M2)[4]*x - coef(M2)[3]), </span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>add=T, lwd=2, col="royalblue")</span> # 2</div>
<div>
<span style="color: lime;">curve(logistic(coef(M2)[4]*x - coef(M2)[3]), add=T, lwd=2, col="darkblue") </span># 3</div>
</div>
<div>
<br />
<br /></div>
<div>
<div>
### 次に、<span style="color: #3d85c6;">推定値の曲線を図示してみる</span>(上のは推定確率でした)</div>
<div>
# 比較対照用にGLMの曲線と比較する</div>
<div>
# 元データ</div>
<div>
<span style="color: lime;">plot(Y1 ~ X, data=D, xlim=c(0,10), ylim=c(0,3))</span></div>
<div>
# GLM</div>
<div>
<span style="color: lime;">curve(3*logistic(coef(M1)[1] + coef(M1)[2]*x), add=T, col="red", lwd=2)</span></div>
<div>
<br /></div>
<div>
# 累積ロジット(全部足し合わせるような累積構造になる)</div>
<div>
<span style="color: lime;">curve(</span><br />
<span style="color: lime;"> 1*(1 - logistic(coef(M2)[4]*x - coef(M2)[1]))</span> # ランク0<br />
<div>
<span style="color: lime;">+ 2*(logistic(coef(M2)[4]*x - coef(M2)[1]) - logistic(coef(M2)[4]*x - coef(M2)[2])) </span># ランク1</div>
<div>
</div>
<div>
<span style="color: lime;">+ 3*(logistic(coef(M2)[4]*x - coef(M2)[2]) - logistic(coef(M2)[4]*x - coef(M2)[3])) </span># ランク2</div>
<div>
</div>
<span style="color: lime;">+ 4*(logistic(coef(M2)[4]*x - coef(M2)[3])) </span># ランク3</div>
<div>
<span style="color: lime;"> - 1, </span># 0, 1, 2, 3なので、1,2,3,4から1を引く</div>
<div>
<span style="color: lime;">add=T, col="blue", lwd=2)</span></div>
</div>
<div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiXCN46_7NhMumYAO_Ol_AOFspUcH3Kl8YJ_qJKU0H8jHVEoUnooucH8f-r0ssI4pwLqL74uciqDNXGIEjAVFbl7g6dR3QdyuS3mKBF1Ik6w-UxwoGpmVJVbdfLDXLICdV9CrFjCa5dNZM/s1600/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588+2014-06-10+18.06.11.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="304" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiXCN46_7NhMumYAO_Ol_AOFspUcH3Kl8YJ_qJKU0H8jHVEoUnooucH8f-r0ssI4pwLqL74uciqDNXGIEjAVFbl7g6dR3QdyuS3mKBF1Ik6w-UxwoGpmVJVbdfLDXLICdV9CrFjCa5dNZM/s1600/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588+2014-06-10+18.06.11.png" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
# 赤:GLM、青:累積ロジット。この単純なケースではほとんど推定値は変わらない。</div>
<br />
<br />
<span style="color: #3d85c6;">#### ついでにベイズ版でも計算</span></div>
<div>
# 下準備、データを累積に変換。1以上、2以上、3以上にまとめる<br />
<span style="color: lime;">rank <- matrix(0, nrow=100, ncol=3) </span># 3はランクの数-1、100はいわゆるN数<br />
<span style="color: lime;">for (r in 1:3) rank[Y1 >= r, r] <- 1 </span><br />
# GLM用の用意していた数値データY1を利用、ランク1~3以上の場合に各列に1を入れる<br />
<br />
<span style="color: lime;">model <- function() {</span><br />
<span style="color: lime;"> for (j in 1:3) { </span># ランクの数-1<br />
<span style="color: lime;"> for (i in 1:100) { </span># いわゆるN数<br />
<span style="color: lime;"> rank[i,j] ~ dbern(p[i,j]) </span># それぞれベルヌーイ分布<br />
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span> logit(p[i,j]) <- alpha[j] + beta*X[i] </span># 切片だけランク毎になってる<br />
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span> }</span><br />
<span style="color: lime;"> alpha[j] ~ dnorm(0, 1E-6) </span># 切片をランク毎に推定<br />
<span style="color: lime;"> }</span><br />
<span style="color: lime;"> beta ~ dnorm(0, 1E-6) </span># 回帰係数は共通<br />
<span style="color: lime;">}</span><br />
<span style="color: lime;"><br /></span>
<span style="color: lime;">data <- list(X=X, rank=rank)</span><br />
<span style="color: lime;">parm <- list(beta=0, alpha=rnorm(3))</span><br />
<span style="color: lime;">source("WBUGS.R") </span># <a href="http://nhkuma.blogspot.jp/p/r2winbugs-macwinewindows-wbugs.html">自作ラッパー関数</a><br />
<span style="color: lime;">out <- wbugs(data, parm, model)</span><br />
<br />
# 3 chains, each with 11000 iterations (first 1000 discarded), n.thin = 10<br />
# n.sims = 3000 iterations saved<br />
# mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff<br />
# beta 1.328 0.166 1.022 1.217 1.318 1.434 1.680 1.001 2200<br />
# alpha[1] -6.377 0.889 -8.217 -6.960 -6.327 -5.770 -4.703 1.002 1100<br />
# alpha[2] -8.645 1.142 -11.040 -9.370 -8.563 -7.851 -6.528 1.002 1800<br />
# alpha[3] -10.518 1.336 -13.331 -11.373 -10.420 -9.585 -8.139 1.001 3000<br />
# deviance 154.787 2.934 151.200 152.600 154.050 156.200 162.000 1.002 1500<br />
# DIC info (using the rule, pD = Dbar-Dhat)<br />
# pD = 4.0 and DIC = 158.8</div>
<div>
<br /></div>
<div>
# BUGS版では、切片が正負が逆になった</div>
<div>
# きちんとコードを書いた結果がこれなので、符号が逆なのはvglmやclmなのだが…じつにややこしい。</div>
Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-39236163392096959232014-06-09T13:28:00.003+09:002014-06-10T18:26:07.473+09:00累積ロジットの汎用Rパッケージ {ordinal}累積ロジット(cumulative logit model)を使う際に、今ひとつ使い勝手のいいRパッケージがないのが気になっていた。<br />
(cf. 累積ロジット:よい、ふつう、わるい、のような段階的な現象についての推定に用いる。等間隔に近ければ二項分布のGLMで構わないだろうけれど、そうでない場合にはこれが適しているようだ)<br />
<br />
例えば、vglmではランダム効果が入れられない。<br />
<br />
mixcatでは、ランダム効果が入れられるが、対数尤度までは出力できても、AICなど情報量規準は出してくれない(手計算すれば良い話ではあるが…)。<br />
<br />
最近、ordinalというパッケージを見つけた。1つのパッケージでGLM版、GLMM版の両方が含まれているし、使用法もglm()やglmer()と同じようにしてくれていて使いやすい。<br />
<br />
ちなみに、GLM版に<span style="color: lime;">stepAIC</span>は使えたが、<span style="color: lime;">dredge</span>はダメだった。<br />
GLMM版では、モデル選択関数は<span style="color: lime;">drop1</span>が使用できた。<br />
<br />
<span style="color: lime;">require(ordinal)</span><br />
<span style="color: lime;">example(ordinal) </span># パッケージで用意されている例を出力<br />
<br />
<span style="color: #3d85c6;"># まずはGLM版</span><br />
ordinl> ## A simple cumulative link model:<br />
ordinl> <span style="color: lime;">fm1 <- clm(rating ~ contact + temp, data=wine)</span><br />
<br />
ordinl> <span style="color: lime;">summary(fm1)</span><br />
formula: rating ~ contact + temp<br />
data: wine<br />
<br />
link threshold nobs logLik <span style="color: #3d85c6;">AIC</span> niter max.grad cond.H<br />
logit flexible 72 -86.49 184.98 6(0) 4.01e-12 2.7e+01<br />
<br />
Coefficients:<br />
Estimate Std. Error z value Pr(>|z|) <span style="color: #3d85c6;"> # 回帰係数部分(この例ではカテゴリカルだが)</span><br />
contactyes 1.5278 0.4766 3.205 0.00135 **<br />
tempwarm 2.5031 0.5287 4.735 2.19e-06 ***<br />
---<br />
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1<br />
<br />
Threshold coefficients: <span style="color: #3d85c6;"> # 各ランクの切片</span><br />
Estimate Std. Error z value<br />
1|2 -1.3444 0.5171 -2.600<br />
2|3 1.2508 0.4379 2.857<br />
3|4 3.4669 0.5978 5.800<br />
4|5 5.0064 0.7309 6.850<br />
<br />
<br />
# 比較対照用にvglmでの結果<br />
<span style="color: lime;">require(VGAM)</span><br />
<span style="color: lime;">vm1 <- vglm(rating ~ contact + temp, family=cumulative(parallel=T), data=wine)</span><br />
<br />
Coefficients:<br />
Estimate Std. Error z value<br />
(Intercept):1 -1.3444 0.50850 -2.6438<br />
(Intercept):2 1.2508 0.43908 2.8487<br />
(Intercept):3 3.4669 0.59711 5.8061<br />
(Intercept):4 5.0064 0.72906 6.8669<br />
contactyes -1.5278 0.47362 -3.2258 # 回帰係数部分の正負がclmと逆なことに注意<br />
tempwarm -2.5031 0.53199 -4.7052<br />
<div>
<br /></div>
<div>
<br />
<span style="color: #3d85c6;"># ordinalに戻って、こちらはGLMM版</span><br />
ordinl> ## A simple cumulative link mixed model:<br />
ordinl> <span style="color: lime;">fmm1 <- clmm(rating ~ contact + temp + (1|judge), data=wine) </span># glmer()同様の構造<br />
<br />
ordinl> <span style="color: lime;">summary(fmm1)</span><br />
Cumulative Link Mixed Model fitted with the Laplace approximation<br />
<br />
formula: rating ~ contact + temp + (1 | judge)<br />
data: wine<br />
<br />
link threshold nobs logLik <span style="color: #3d85c6;">AIC</span> niter max.grad cond.H<span style="color: #3d85c6;"> # ちゃんとAICも算出する</span><br />
logit flexible 72 -81.57 177.13 331(996) 1.05e-05 2.8e+01<br />
<br />
Random effects: <span style="color: #3d85c6;"># ランダム効果</span><br />
Groups Name Variance Std.Dev.<br />
judge (Intercept) 1.279 1.131 <br />
Number of groups: judge 9<br />
<br />
Coefficients:<br />
Estimate Std. Error z value Pr(>|z|) <br />
contactyes 1.8349 0.5125 3.580 0.000344 ***<br />
tempwarm 3.0630 0.5954 5.145 2.68e-07 ***<br />
---<br />
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1<br />
<br />
Threshold coefficients:<br />
Estimate Std. Error z value<br />
1|2 -1.6237 0.6824 -2.379<br />
2|3 1.5134 0.6038 2.507<br />
3|4 4.2285 0.8090 5.227<br />
4|5 6.0888 0.9725 6.261<br />
<div>
<br />
(<a href="http://nhkuma.blogspot.jp/2014/06/glm.html">次の記事に続く</a>)</div>
</div>
Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-90497420818566528872014-03-15T13:19:00.001+09:002014-03-15T13:19:39.225+09:00生態学会2014広島大会で発表します<span style="background-color: rgba(255, 255, 255, 0);">広島の生態学会に来ています。発表のネタの準備中に大きなミスを見つけて取り返したりで時間に追い詰められグッタリしてますorz</span><br />
<span style="background-color: rgba(255, 255, 255, 0);">発表は明後日17日のポスター、国内温帯域の藻場の衰退とサンゴの分布拡大の話です。温暖化影響や生物多様性データベースのGIS化などに興味のある方には興味の持てる内容かと思います。来聴歓迎します。</span><br />
<span style="background-color: rgba(255, 255, 255, 0);"><a href="http://www.esj.ne.jp/meeting/abst/61/PB3-083.html">http://www.esj.ne.jp/meeting/abst/61/PB3-083.html</a></span>Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-58402067586946583362014-02-20T15:36:00.001+09:002014-03-25T10:44:49.896+09:00(旧サイトより移行)過分散データ:GLM負の二項分布、GLMMによる推定をAICでモデル選択することは可能か?# 2014.02.20追記:下記のlmer()は現在はglmer()に相当します<br />
<br />
かつて、数値実験をしてみて「できないようだ」と書いたのですが、<span style="color: orange;">GLMM側(lmerやglmmML)のAIC計算方法をGLM側(glmやglm.nb)へと合わせるように修正すれば可能</span>であることを確かめました。じつは、両者でAICの計算方法が異なっていたことが分かりました。道理でいつもGLMMのAICが小さ過ぎるわけです…。<br />
<br />
★以前アップしていた数値実験を大幅に修正・加筆しておきました(cf.「<a href="http://nhkuma.blogspot.jp/2014/02/glmglmm.html">過分散データ:GLM負の二項分布、GLMMによる解析の比較</a>」)<br />
<br />
なお、glmやglm.nb と、lmerやglmmMLのAICの計算方法の違いを検証するところまでは、検索すると複数のブログで見つかります。<br />
<br />
ところがさらに検証してみたところ、<span style="color: orange;">逆にglmやglm.nbのAICの方をlmerやglmmMLのAIC計算方法へと修正すると、不当にglm、とくにglm.nbのAICが小さくなり、誤ったモデル選択となる確率が大幅に増大することも分かりました(ほぼ逆転してしまう!)</span>。もっとも、ちゃんと数値実験をやっているとはいえ、なぜそうなるのか理論的裏付けがよく分からないので、なにか思いついたら検証したいと思います。<br />
<br />
ちなみに、まったくのポアソン分布データのパラメータ推定でも、glm、glm.nb、lmerはいずれも真の値をほぼ推定できていました。過分散問題がけっこう厄介なことを思うと、もはや glm(, family=poisson) を紹介する意義はほとんど無いのかもしれません…。<br />
<br />
<br />
<br />
数値実験では、乱数発生させた、<span style="color: #3d85c6;">ポアソン分布データと、負の二項分布データ、ポアソン分布に正規乱数ノイズを加えたデータの三者で、GLM(ポアソン分布:glm)とGLM(負の二項分布:glm.nb)、GLMM(ポアソン分布:lmer)が真の値を推定できているかをクロスチェックし、同時にAICが正しい推定をちゃんと反映できているかをもチェック</span>しました。lmerの方は本来のGLMMの使い方とは異なり、各データを異なるid(1データ1id)としているので過分散の対処に用途は限定しています(こういう使い方は、help(lmer)や、Warton & Hui (2011) Ecology 92:3–10、和文では粕谷先生のGLM本「一般化線形モデル」で紹介されています)。<br />
<br />
混乱しやすいのでメモ的に整理します。<br />
まず、AIC()を用いた時の、glmやglm.nbと、lmerやglmmMLの計算方法の違い:<br />
glmやglm.nbでは、<br />
AIC = deviance + 2*推定するパラメータ数<br />
lmerやglmmMLでは、<br />
AIC = residual deviance + 2*推定するパラメータ数<br />
<br />
また、deviance()は、いずれの場合もresidual devianceを計算する。<br />
一方、logLik()では、glmやglm.nbではdeviance/(-2)、lmerやglmmMLではresidual deviance/(-2)<br />
<br />
<br />
lmerやglmmMLのAICをglmの方へ合わせるには、下記のようにする(上述の理由により、逆はやめておきましょう)。<br />
<br />
# 例えば、こんな2つのモデルがあるとき(id=1~N数)<br />
<span style="color: lime;">pois <- glm(Y ~ X, family=poisson)</span><br />
<span style="color: lime;">pois.m <- lmer(Y ~ X + (1|id), family=poisson)</span><br />
<br />
glmのAICはふつうに、AIC(pois)、で求められる。<br />
<span style="color: #3d85c6;">GLMMのAICは、上述の違いに基づくと、こうすればglmのと同じように求められる:</span><br />
<span style="color: lime;">AIC(pois.m) - deviance(pois) -2*logLik(pois) </span><br />
# (residual deviance + 2*パラメータ数) - residual deviance + deviance = deviance + 2*パラメータ数<br />
#(ランダム効果を抜く以外は同等であるglmモデルを用意してやる必要がある)<br />
<br />
# ちなみに、<span style="color: #3d85c6;">AICがそのように変換できることの検証として、GLMMでidを全部 1 にしてやると(1グループのみ)、グループ間の分散は 0 と推定されます。切片と回帰係数の推定も、単なるGLMのそれと同一の推定結果となります。ただし、分散値は推定するパラメータ数としてカウントされるので、変換後のGLMMのAICはGLMのAICよりもキッカリ2だけ増加した値となります。</span>これは、AICの計算式にある、2*パラメータ数(ここでは、2*1)、の部分に相当します。<br />
<br />
しかし、なぜGLMMのパッケージがresidual devianceに基づいたAIC計算にしているのか、その意図が何なのかわかりません。こういう計算方法の修正をしないでくれと言っているようにも見えます。さしあたり、過分散データへの対処という限定的な使用法について数値実験をしてみた限りでは問題は無いようです。<br />
<br />
こういう煩わしさを思うと、もはや個人的にはBUGSに全面移行してしまおうかなどと考えたりします…Stanの発展も待ち遠しいです。<br />
<br />
<br />
# (当時いただいたコメント)<br />
2012/11/15(木) 23:39:11 | URL | 高橋<br />
興味深く読ませていただきました。<br />
ちょっと気になったのですが、<br />
># (residual deviance + 2*パラメータ数) - residual deviance + deviance = deviance + 2*パラメータ数<br />
この式の2つのresidual devianceは別物だと思います。<br />
residual deviance = (当該モデルのdeviance) - (飽和モデルのdeviance)<br />
なので、<br />
>- deviance(pois) -2*logLik(pois)<br />
は、<br />
- ((ポアソンモデルのdeviance) - (ポアソンモデルの飽和モデルのdeviance)) + (ポアソンモデルのdeviance)<br />
となります。<br />
したがって、<br />
>AIC(pois.m) - deviance(pois) -2*logLik(pois)<br />
は、<br />
(混合モデルのdeviance) - (混合モデルの飽和モデルのdeviance) + 2*パラメータ数 + (ポアソンモデルの飽和モデルのdeviance)<br />
となり、「混合モデルの飽和モデル」が「ポアソンモデルの飽和モデル」と等しい場合(多分そうなのだと思いますが)、混合モデルのdevianceに基づいてAICを求めているということになると思います。<br />
glm.nb()の関数定義を見ると、推定されたThetaのもとで飽和モデルのdevianceを計算しているようです。この値は、ポアソンモデルにおける飽和モデルのdevianceよりずっと大きくなるので、residual devianceを用いてAICを計算すると、負の二項分布の常勝となるのでしょう。<br />
ちなみに飽和モデルの対数尤度は、応答変数をyとすれば、<br />
sum(dpois(y, y, log = TRUE))<br />
で求められます。負の二項分布の場合は、glm.nb()が返すThetaパラメータを使って、<br />
sum(dnbinom(y, mu = y, size = Theta, log = TRUE))<br />
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<div>
> 高橋さま </div>
<div>
コメントありがとうございます、おかげで理解が深まりました。この検証に興味を持っていただき嬉しいです。 </div>
<div>
<br /></div>
<div>
ポアソンモデルと、ポアソン混合モデルの飽和モデルが等しい時にのみ成り立ちうるというご指摘、たしかにその通りでした。 </div>
<div>
<br /></div>
<div>
これはもう数値実験を見る限りそうなるとしか答えようがないのですが、1グループのみの混合モデルの変換後AICから2を引いた値とポアソンモデルのAICが一致する、というのが一応の証拠になっていると考えています(こちらの検証コード掲載は省いてしまっていますが)。 </div>
<div>
<br /></div>
<div>
glm.nbの飽和モデルの求め方がそうなっているとは気づいていませんでした。言われてみて改めて関数定義をチェックするとたしかにそれらしきコードを見つけることができました。 </div>
<div>
<br /></div>
<div>
混合モデルならば、グループを1つだけにするなどにより、ポアソンモデルとのすり合わせを試みることができましたが、負の二項モデルで同じようなすり合わせをするとなると…glm関数でThetaを1に固定してみるか…また追々試してみたいと思います。</div>
</div>
<div>
<br /></div>
Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-18857804757928159962014-02-20T15:30:00.000+09:002014-02-20T15:30:03.563+09:00(旧サイトより移行)連続量を単位あたりに直してGLMする場合の対処:gaussian("log")、Gamma("log")、offsetなどなど)# 2014.02.20 旧サイトを閉じるため、移植しました。<br />
# なお、下記のようなケースではゴンペルツ曲線などを用いた成長モデルを使った方がいいかなと最近思ってます。<br />
<br />
最近、実験系の統計を引き受けていて気になったので、連続量を単位あたりでGLM推定する方法について検証してみました。<br />
<br />
単位あたりの量をGLMで解析する時、個体数や頻度などの整数の場合の対処はoffset項の利用で解決するのが常套手段ですね(密度や○×率のような割り算した値ではなく、元の値はいじらずに単位量を係数1として説明変数に加える。当ブログでもかつて紹介:http://nhkuma269.blog77.fc2.com/blog-entry-9.html)。<br />
例えば、同じく個体密度0.1でも、1/10と100/1000とでは意味合いが違うが、割り算すると両者は同一密度として扱われてしまう。前者では1個体の増減が大きな誤差を生むが、後者ではほとんど影響なし。<br />
<br />
<span style="color: #3d85c6;">では、元々が連続量のものを単位あたりの量にする場合も同様にoffset項による対処がよいのだろうか?例えば、サンゴの枝あたりのクロロフィル量とか、成長量(実験後重量 / 実験前重量)のような量は、必ず単位量(この場合は枝)のバラツキによる誤差が出てしまう(なるべく条件は揃えるように努力はしているだろうが)。これを割り算してしまうと誤差が統計結果に悪影響するだろう。</span><br />
<br />
CrawleyのR統計本によると、単位量は共変量として説明変数に加えるか、log(目的の量 / 単位量)というように、割り算した上で対数変換せよ。と書かれていた。前者は要因の1つとしてカウントするならばそれでいいとして、後者はすでに古典的な対処法でしょう…(本自体がすでに一時代前のものなので仕方ないですが)。<br />
<br />
ついでに、論文などで時々見掛ける、family=gaussian(link="log")、正規分布でリンク関数が対数という一見すると対数正規分布っぽいけれど違うらしい(対数正規分布の分散は平均と共に増大するようだ→平均や分散の式をチェックされたし)、こいつの性質もチェックしてみます。もう一つおまけに、古典統計で対数変換する際によく使うlog(x + 1)変換もチェックしてみました。ただし、これはlog(0)が計算できない問題を回避する目的に限定します。<br />
<br />
<br />
<span style="color: #3d85c6;">比較するモデルは下記の6バージョン。これを応答変数が対数正規分布とガンマ分布の場合のそれぞれについて一次回帰式のパラメータ推定(切片と回帰係数)をします:</span><br />
(色は下記の図とおおよそ対応)<br />
<span style="color: #cccccc;">m11:正規分布モデル、割合にした上でlog(x + 1)の対数変換という古典統計の常套手段(グラフ上では黒ライン)</span><br />
<span style="color: red;">m12:正規分布モデル、割合にした上でlog(x + 1e-10)変換、1の代わりにごく小さい小数(1.0 x 10^(-10))でlog(0)を回避</span><br />
<span style="color: lime;">m13:正規分布モデル、割合にしているが、変数変換はせず、連結関数を対数に指定(以下、log(0)回避には1e-10を使用)</span><br />
<span style="color: blue;">m14:正規分布モデル、単位量をoffset(log(off))として係数1の説明変数に加え、連結関数を対数に指定</span><br />
<span style="color: cyan;">m23:ガンマ分布モデル、割合にしているが、変数変換はせず、連結関数を対数に指定(ガンマ分布は正の値しか取れないので、1と2の変数変換のモデルは作れない)</span><br />
<span style="color: magenta;">m24:ガンマ分布モデル、単位量をoffset(log(off))として係数1の説明変数に加え、連結関数を対数に指定</span><br />
<br />
<br />
結論から言うと、ケースバイケースな複雑な結果になりました…。<span style="color: orange;">実際の利用では、m13とm23、またはm14とm24をAICで比較するのがよいでしょう。</span><br />
<span style="color: orange;">・やはりデータの対数変換は止めた方がよい、連結関数:logを用いるべき</span>です。対数変換したモデルでは推定が大きく乱れる場合があり信頼できない(図1_1, 1_2のm11黒、m12赤)。<br />
・データが対数正規分布の時、割合の正規分布モデル(m13緑)は、offsetの正規分布モデル(m14青)よりも推定の分散が小さかった(図1_1、2_1)。ただし、データがガンマ分布の時は、m14の方がm13よりも推定の分散が小さかった(図1_2、2_2)。両者の違いはさほど大きくなかった。<br />
・ガンマ分布モデルでは、割合(m23水色)と、offset(m24ピンク)とは推定値がまったく同じ!(そのため、m23の水色は完全にマスクされている)<br />
・つまり、この数値実験の限りでは、予想に反して<span style="color: orange;">連続量を割合や比率にした量をGLMで解析する際、割り算をした値を用いても大した問題はない</span>ことになる。offsetの使用による論文中での説明のややこしさや、記述可能なモデルの可塑性などを考えると…割り算していい気がしてきました。しかし、なぜ割り算にしても大丈夫だったのか、背景にある数学的なロジックは朧気なままです…たぶん、元の値と割算値とで確率分布が変わらないからだと思います(よくある整数値の割り算の場合は、小数点を取るようになるので本来のポアソンや二項分布が使用不可になるが)。例外があるとすれば、単位量の方にも確率誤差(その値が属している確率分布からのズレ)がある場合でしょう。その場合はoffsetすべきということかと。<br />
・ちなみに、<span style="color: orange;">常套手段なlog(x + 1)変換(m11黒)の推定は危険です。とくに平均値の小さい推定では危険(図1_1、1_2)。log(0)の回避には1の代わりにごく小さな値(1e-10くらい)を足すのがよいでしょう。</span>ただし、古典統計で整数を対数変換する場合は、平均と分散の関係を調整する意図があるので + 1 のままで)。<br />
<br />
なお、ここでやっているのは入り口と出口がちゃんと一致するかを確認するだけの数値実験に過ぎません。参考にする際には自己責任でお願いします。<br />
<br />
<br />
<br />
<span style="color: #3d85c6;"># 以下、こんな関数でGLM推定を1000回繰り返し、パラメータ推定の精度をチェックしました。</span><br />
<span style="color: #3d85c6;"># 平均 exp(alpha + beta*X)という一次回帰を考え、こちらで指定したalphaとbetaを用いて、対数正規分布とガンマ分布のデータを発生させ、各モデルによってalphaとbetaを逆推定し、初めに指定した値を再現できるかどうかチェックします。</span><br />
# alpha, beta, distributionは1が対数正規分布、2がガンマ分布、n.itrは繰り返し数を表します。<br />
<br />
cont.skew <- function(alpha, beta, distribution, n.itr) {<br />
estim <- numeric(0) # 後でデータをくっつけるためのイントロンみたいなもの<br />
set.seed(1)<br />
X <- rep(c(1:10)*0.1, each=10)<br />
off <- runif(length(X), min=1, max=10) # 下記のoffsetで用います。数字に特に意味はなし<br />
mean <- exp(alpha + beta*X)*off<span class="Apple-tab-span" style="white-space: pre;"> </span># 平均:log(Y/off) ~ alpha + beta*X、なので。<br />
sd <- exp(0.5)<br />
for (n in 1:n.itr) {<br />
Y <- rbind(rlnorm(100, meanlog=log(mean), sdlog=log(sd)), # 対数正規分布<br />
rgamma(100, shape = mean^2/sd^2, scale = sd^2/mean)) # ガンマ分布(cf. rgamma()のhelp)<br />
Y <- Y[distribution, ]<br />
D <- data.frame(X, Y, off)<br />
<span style="color: #cccccc;">m11 <- glm(log(Y/off + 1) ~ X, family=gaussian(link="identity"), D) # 正規分布モデル(割合、対数変換1)</span><br />
<span style="color: red;">m12 <- glm(log(Y/off + 1e-10) ~ X, family=gaussian(link="identity"), D) # 正規分布モデル(割合、対数変換2)</span><br />
<span style="color: lime;">m13 <- glm((Y/off + 1e-10) ~ X, family=gaussian(link="log"), D) # 正規分布モデル(割合、連結関数:対数)</span><br />
<span style="color: blue;">m14 <- glm((Y + 1e-10) ~ X + offset(log(off)), family=gaussian(link="log"), D) # 正規分布モデル(オフセット、連結関数:対数)</span><br />
<span style="color: cyan;">m23 <- glm((Y/off + 1e-10) ~ X, family=Gamma(link="log"), D) # ガンマ分布モデル(割合、連結関数:対数)</span><br />
<span style="color: magenta;">m24 <- glm((Y + 1e-10) ~ X + offset(log(off)), family=Gamma(link="log"), D) # ガンマ分布モデル(オフセット、連結関数:対数)</span><br />
estim2 <- c(rbind(coef(m11), coef(m12), coef(m13), coef(m14), coef(m23), coef(m24)))<br />
estim <- rbind(estim, estim2)<br />
} # forループ、ここまで<br />
par(mfcol=c(1,2))<br />
plot(density(estim[,1]), lwd=4, xlim=c(alpha-0.5, alpha+2), ylim=c(0,3), main="alpha (intercept)")<br />
for (i in 2:6) { lines(density(estim[,i]), lwd=4, col=i)<br />
abline(v=alpha) } # alpha値の図示<br />
plot(density(estim[,7]), lwd=4, xlim=c(beta-1, beta+1), ylim=c(0,3), main="beta (coefficient)")<br />
for (j in 2:6) { lines(density(estim[,j+6]), lwd=4, col=j)<br />
abline(v=beta) } # beta値の図示<br />
estim.m <- matrix(apply(estim, 2, mean), ncol=2) # alpha, betaの推定値を計算<br />
colnames(estim.m) <- c("alpha", "beta")<br />
estim.m} # 推定値を表示(cont.skew関数、ここまで)<br />
<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgRknElNXRHy7YfsJa6RfFMow96gkLhYm-D1G4g8OR4gq5aELQDfZ-P8xaMhGOZ2RNXikXKv-JWUP0pVg_qfTyIoWnfzS5Eu7_0D1yS8zktKhnP61LRZQNIYenCg8OvgAFkj8_w4sSWEng/s1600/20120529161852dd0.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgRknElNXRHy7YfsJa6RfFMow96gkLhYm-D1G4g8OR4gq5aELQDfZ-P8xaMhGOZ2RNXikXKv-JWUP0pVg_qfTyIoWnfzS5Eu7_0D1yS8zktKhnP61LRZQNIYenCg8OvgAFkj8_w4sSWEng/s1600/20120529161852dd0.jpg" height="207" width="320" /></a></div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<span style="color: #3d85c6;"># 図1_1: 対数正規分布のパラメータ推定、平均値が小さい場合</span><br />
cont.skew(alpha=-1, beta=1, n.itr=1000, distribution=1)<br />
# m11:黒、m12:赤、m13:緑、m14:青、m23:水色、m24:ピンク(水色はピンクと完全に重なって表示されない)<br />
# m12赤、m23&m24ピンク、m13緑、の順に推定がよかった。m11黒は大きくずれた。<br />
# 推定値(1~6の順にm11~m24)、m23とm24(一番下の2つ)は完全に推定値が一致!<br />
# alpha beta<br />
# [1,] 0.3131626 0.3959152<br />
# [2,] -1.0009638 0.9999089<br />
# [3,] -0.8811454 1.0011294<br />
# [4,] -0.8879904 1.0073653<br />
# [5,] -0.8793282 0.9996438<br />
# [6,] -0.8793282 0.9996438<br />
<br />
<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiSt93oLpxBPm28qgvo9qqliZNs7KIyjDLgwpBSn_pOOJyoS4x5G5LMVni7xqd7Sv1ar5rGKRFUEnM9GniSDG45ZKN4LovwoEPbI4VkCooFCxmlXo86Z_CBllxB_68dY9wHY6KxzFI-ems/s1600/2012052916185183c.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiSt93oLpxBPm28qgvo9qqliZNs7KIyjDLgwpBSn_pOOJyoS4x5G5LMVni7xqd7Sv1ar5rGKRFUEnM9GniSDG45ZKN4LovwoEPbI4VkCooFCxmlXo86Z_CBllxB_68dY9wHY6KxzFI-ems/s1600/2012052916185183c.jpg" height="207" width="320" /></a></div>
<br />
<span style="color: #3d85c6;"># 図1_2:ガンマ分布のパラメータ推定、平均値が小さい場合</span><br />
cont.skew(alpha=-1, beta=1, n.itr=1000, distribution=2)<br />
# m11:黒、m12:赤、m13:緑、m14:青、m23:水色、m24:ピンク(水色はピンクと完全に重なって表示されない)<br />
# m14青、m13緑、m23&m24ピンク、の順に推定がよかった。m11黒とm12赤は大きくずれた。<br />
# なぜかガンマ分布の推定なのにガンマ分布モデル(m23&m24)が最強にならない…<br />
<br />
<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhSUaueEooiMAbLWiL5BzbB86D8tOuzjGTQS0oRMnz9l7Hnx9yDonzvRx7iUbi5GEL4Fev1fTTNnu_9F2PEQR3cOw1-HYgGOocqkj76sU0JbNLXWDJ-5Wia4nRxSIsUW-78UMU7hxJaNoE/s1600/20120529161850869.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhSUaueEooiMAbLWiL5BzbB86D8tOuzjGTQS0oRMnz9l7Hnx9yDonzvRx7iUbi5GEL4Fev1fTTNnu_9F2PEQR3cOw1-HYgGOocqkj76sU0JbNLXWDJ-5Wia4nRxSIsUW-78UMU7hxJaNoE/s1600/20120529161850869.jpg" height="207" width="320" /></a></div>
<br />
<span style="color: #3d85c6;"># 図2_1:対数正規分布のパラメータ推定、平均値が大きい場合</span><br />
cont.skew(alpha=2, beta=1, n.itr=1000, distribution=1)<br />
# m11:黒、m12:赤、m13:緑、m14:青、m23:水色、m24:ピンク(水色はピンクと完全に重なって表示されない)<br />
# 平均値が大きくなると、モデルによる推定の違いは小さくなる(中心極限定理による正規分布への収束)<br />
# m12赤、m23&m24ピンク、m11黒、の順に推定がよかった。<br />
<br />
<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiBpRpAUoPq2L126vHoTgzgz84U8dz7ShIG5ip7aDVbDDwDwtJtAyyD-LbR4Ltb9948v_yfx3pvEZK6GUxTGXa9wua00TCqdmULQfmL5eEZWFjmbfc6kxqPAx9d2o4keeu_eW62Y6kqy20/s1600/20120529161849bc9.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiBpRpAUoPq2L126vHoTgzgz84U8dz7ShIG5ip7aDVbDDwDwtJtAyyD-LbR4Ltb9948v_yfx3pvEZK6GUxTGXa9wua00TCqdmULQfmL5eEZWFjmbfc6kxqPAx9d2o4keeu_eW62Y6kqy20/s1600/20120529161849bc9.jpg" height="207" width="320" /></a></div>
<br />
<span style="color: #3d85c6;"># 図2_2:ガンマ分布のパラメータ推定、平均値が大きい場合</span><br />
cont.skew(alpha=2, beta=1, n.itr=1000, distribution=2)<br />
# m11:黒、m12:赤、m13:緑、m14:青、m23:水色、m24:ピンク(水色はピンクと完全に重なって表示されない)<br />
# 平均値が大きくなると、モデルによる推定の違いは小さくなる(中心極限定理による正規分布への収束)<br />
# m11黒以外は真の値にほぼ収束Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com2tag:blogger.com,1999:blog-6412094354437530509.post-50644067364501745312014-02-20T15:15:00.002+09:002014-06-24T20:20:58.673+09:00(旧サイトより移行)過分散データ:GLM負の二項分布、GLMMによる解析の比較## 2014.02.20 追記:旧サイトを閉じるため、このページを移植しました。<br />
なお、現在は lmer()は正規分布専用になっており、その他の分布ではglmer()関数を使用します。<br />
<br />
##(2012.10.17 追記:GLM関数群とGLMM関数群との間でのモデル選択について、数値実験結果を追加・修正しました)<br />
<br />
# 生態学会で話していて、<span style="color: #3d85c6;">過分散データの取り扱いが話題に上った。負の二項分布を使えという意見と、GLMMを使えという意見と両方があるが、実際どっちの方がよいというのはあるのだろうか?</span>(cf. 2014.6.24追記: <a href="http://nhkuma.blogspot.jp/2014/02/r-winbugs.html">GLM負の二項分布とは実態はランダム効果がガンマ分布となったGLMM</a>)<br />
<br />
<span style="color: #3d85c6;"># 過分散 overdispersion:ポアソン分布や二項分布のモデルに当てはめる場合、これらの分布型が平均と分散が表裏一体なのだが(ポアソンでは平均=分散)、実際のデータは大抵はそうでないので分散が過剰となることが多い、という問題。</span><br />
# GLMのsummary出力結果で、下の方に、こんなのがあって…<br />
# Residual deviance: 5094.7 on 71 degrees of freedom<br />
# …ここで、5094.7/71の値が、1.5を超えているならば、何らかの対処をした方がよいだろう(cf. モデル選択ペンギン本、Zuur et al (2009) Mixed effects models and extensions in Ecology with R. Springer。これのp. 225-6には"1.5"が目安と書かれている)<br />
<br />
<span style="color: #3d85c6;"># 過分散しているデータとは下記のように裾野が片方に広く延びている(上:ふつうのポアソン分布、中:負の二項分布、下:正規ノイズを加えたポアソン分布)。</span>負の二項分布はゼロのあたりに集中、正規ノイズは正負にばらつく。<br />
<br />
par(mfcol=c(3,1)) # グラフ領域を縦に3分割<br />
hist(rpois(10000, lambda=exp(1)), freq=F, xlim=c(0,100), breaks=seq(0,20,by=2), col=1) # ポアソン分布(平均 exp(1))<br />
hist(rnegbin(10000, mu=exp(1), theta=1/4), freq=F, xlim=c(0,100), breaks=seq(0,10000,by=2), col=2) # 負の二項分布(平均 exp(1)、θ 1/4。θはガンマ分布のshapeパラメータに相当、shape=平均^2 / 分散)<br />
hist(rpois(10000, lambda=exp(1 + rnorm(10000, mean=0, sd=4^0.5))), freq=F, xlim=c(0,100), breaks=seq(0,10000,by=2), col=3) # ポアソン分布(平均 exp(1))+正規ノイズ(平均0、分散4)<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjnaGfisJd8WBmqVKi2v0Z1vKEctA0RzwvQ0205G9l3hO80N9iCLz7yc7bgZ9Mti5w-mfJCZWKxjYSvLlPlw9I0p4CLHRXCBL0Dq6BNj9j1rnm1hR8n8Thb7A7L7ujDwv4QLqsVazdUpkA/s1600/2012061916044117a.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjnaGfisJd8WBmqVKi2v0Z1vKEctA0RzwvQ0205G9l3hO80N9iCLz7yc7bgZ9Mti5w-mfJCZWKxjYSvLlPlw9I0p4CLHRXCBL0Dq6BNj9j1rnm1hR8n8Thb7A7L7ujDwv4QLqsVazdUpkA/s1600/2012061916044117a.jpg" height="320" width="176" /></a></div>
<br />
<br />
<br />
<br />
<span style="color: #3d85c6;"># では、ポアソン分布のデータと、負の二項分布をしたデータ、ポアソン分布に正規乱数ノイズを加えたデータを人工的に作成し、それぞれGLM (poisson)、GLM(負の二項分布)、GLMM(データ数分の個体差があると見なす)で推定した結果を比べてみよう(実行するためのプログラムは末尾に)。</span><br />
<br />
<br />
OD(alpha=2, beta=0.5, siml=100, distribution=1)<span style="color: #3d85c6;"> # データがポアソン分布の時</span><br />
#(マシンパワーによってはけっこう時間が掛かる、アラートも沢山出るが無視)<br />
# ランダム効果がN数と同じ数だけあるよ!?というアラートが大量に出てくるが、推定はちゃんとできているので無視<br />
# Number of levels of a grouping factor for the random effects is *equal* to n, the number of observations<br />
<br />
<span style="color: #3d85c6;"># 黒:GLM (poisson)、赤:GLM(負の二項分布)、緑:GLMM。100回分の推定結果</span><br />
#(3本のラインがあるはずですが完全に重なっていて、視覚的には識別不能)<br />
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjains2CIBAS907Z_Lp8_Sm1RKUlrggI4l2V35VhNST-okeNoDXSc2V3Bk7AXmgEP8ybHGoV9crAFXuV0fe91Fyfqj8U7nrdpEBrA95YFRVG8TEuR2Fbg_6ldwoLdwMlUaQX6j877_Wf5c/s1600/pois.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjains2CIBAS907Z_Lp8_Sm1RKUlrggI4l2V35VhNST-okeNoDXSc2V3Bk7AXmgEP8ybHGoV9crAFXuV0fe91Fyfqj8U7nrdpEBrA95YFRVG8TEuR2Fbg_6ldwoLdwMlUaQX6j877_Wf5c/s1600/pois.jpg" /></a></div>
<br />
# models: glm.pois glm.nb glmm.pois<br />
# aic.d 710.5256 712.1975 712.2007 # devianceに基づくAIC<br />
# prob.d 0.9500 0.0400 0.0100 # 各モデルが選択される確率<br />
# aic.rd 102.8072 101.6366 104.4823 # residual devianceに基づく AIC<br />
# prob.rd 0.6000 0.4000 0.0000 # 各モデルが選択される確率<br />
<br />
<span style="color: #3d85c6;"># GLM (poisson)、GLM(負の二項分布)、GLMMともalpha、betaの推定値はほぼ真の値となった(ということは、単なるGLM(poisson)を使う理由はないのではないか!?)。AIC(deviance)によるモデル選択は確率0.95で真のモデルであるGLM(ポアソン分布)を選択したが、AIC(residual deviance)ではAICの推定値が僅差で負の二項分布モデルが最小だし、誤って負の二項分布モデルを選択する確率が0.4もあった。とはいえ、いずれのモデルもほぼ真の値を推定できているので問題なしか!?</span><br />
<br />
<br />
<br />
OD(alpha=2, beta=0.5, siml=100, distribution=2) <span style="color: #3d85c6;"># データが負の二項分布の時</span><br />
<span style="color: #3d85c6;"># 黒:GLM (poisson)、赤:GLM(負の二項分布)、緑:GLMM。100回分の推定結果</span><br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjgwXl5AVRvQvhceWcbHoQnX151v3CBeJrvaW3lSFeR_mNLtBwhnqfnHAmGHcprzNR1jGoesNhRJ6Q5Bs-s_g1L3V_xl4_0lb5RNOhGLMIoAzLx7fvZ1bzYDPJjXC2jo3IUvkW55qCWsMI/s1600/negbin.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjgwXl5AVRvQvhceWcbHoQnX151v3CBeJrvaW3lSFeR_mNLtBwhnqfnHAmGHcprzNR1jGoesNhRJ6Q5Bs-s_g1L3V_xl4_0lb5RNOhGLMIoAzLx7fvZ1bzYDPJjXC2jo3IUvkW55qCWsMI/s1600/negbin.jpg" height="236" width="320" /></a></div>
<br />
# models: glm.pois glm.nb glmm.pois<br />
# aic.d 36960.67 902.9074 922.3840 # devianceに基づく AIC<br />
# prob.d 0.00 1.0000 0.0000 # 各モデルが選択される確率<br />
# aic.rd 36572.47 121.0299 534.1787 # residual devianceに基づく AIC<br />
# prob.rd 0.00 1.0000 0.0000 # 各モデルが選択される確率<br />
<br />
<span style="color: #3d85c6;"># GLM (poisson)は両方向に大きく広がった。GLM(負の二項分布)はalphaもbetaも真の値に近かった。GLMMはbetaの値は真の値に近かったが分布の裾が広かった、alphaは盛大にズレた。</span><br />
<span style="color: #3d85c6;">モデル選択はすべてGLM(負の二項分布)となった。</span><br />
<span style="color: #3d85c6;"><br /></span>
<span style="color: #3d85c6;"><br /></span>
<span style="color: #3d85c6;"><br /></span>
<span style="color: #3d85c6;"># ポアソン分布に正規乱数ノイズを加えたデータを同様に推定してみる</span><br />
OD(alpha=2, beta=0.5, siml=100, distribution=2)<br />
<br />
<span style="color: #3d85c6;"># 黒:GLM (poisson)、赤:GLM(負の二項分布)、緑:GLMM。100回分の推定結果</span><br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjQFONeV4JFt5TmZ5hgn_l_A0-y7YzCquQx5p-z_mu2uv2mQCEgfyg18-j9d2WPGrFlVcNznisSg33vbhGcoVlwF5EJLDcje9wx641xFoooJQEvgnrN9N_2LXa8J3-aparrdLBuiAAak-o/s1600/poisnorm.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjQFONeV4JFt5TmZ5hgn_l_A0-y7YzCquQx5p-z_mu2uv2mQCEgfyg18-j9d2WPGrFlVcNznisSg33vbhGcoVlwF5EJLDcje9wx641xFoooJQEvgnrN9N_2LXa8J3-aparrdLBuiAAak-o/s1600/poisnorm.jpg" height="236" width="320" /></a></div>
<br />
<br />
# models: glm.pois glm.nb glmm.pois<br />
# aic.d 845220.9 1309.4785 1278.7899 # devianceに基づく AIC<br />
# prob.d 0.0 0.0100 0.9900 # 各モデルが選択される確率<br />
# aic.rd 844616.2 136.0312 674.1263 # residual devianceに基づく AIC<br />
# prob.rd 0.0 1.0000 0.0000 # 各モデルが選択される確率<br />
<br />
<span style="color: #3d85c6;"># GLM (poisson)は真の値から大きく外れた。GLM(負の二項分布)はalphaは盛大にズレた、betaは真の値に近かったが裾野が広かった。GLMMはalphaもbetaも真の値に近かった。AIC(deviance)によるモデル選択は0.99の確率で真の母集団のモデルGLMMを選択したが、AIC(residual deviance)によるモデル選択はすべてのケースで誤ってGLM(負の二項分布)を選択してしまった…。</span><br />
<span style="color: #3d85c6;"><br /></span>
<span style="color: #3d85c6;"><br /></span>
<span style="color: #3d85c6;"><br /></span>
<span style="color: #3d85c6;"># 結論として、過分散データと言っても一辺倒ではなく、負の二項分布、正規ノイズ、どちらに近いかによって結果が変わってくるという、ごく当然の結果となった。glm.nb、GLMMのどちらがよいかは、結局データ依存ということになりました。</span><br />
<br />
<span style="color: #3d85c6;"># 選択基準としては…glm.nb、GLMMでAICを基にモデル選択は、一手間加えれば可能であるようだ。また、glmとglm.nbの方をresidual devianceに基づくAICに修正する計算も試してみたところ(2*logLik(pois) + deviance(pois) + AIC(pois))、誤ってglm.nbモデルが選択される確率が大きく増大することがわかった。なので、</span><span style="color: #e69138;">GLMとGLMMのモデルをAICで比較する場合には、GLMMのAICをdevianceベースのAICに修正して用いる必要がある</span><span style="color: #3d85c6;">(-2*logLik(pois) - deviance(pois) + AIC(pois.mm))ということですね。</span>(注:poisはglm(,family=poisson)のモデル、pois.mmはlmer(,family=poisson)のモデル、両者の違いはランダム効果の有無のみ → 下記のプログラム中にモデル式あり)<br />
<br />
<br />
# しかし、lmerやglmmMLのデフォルトがresidual devianceに基づくAICを計算していること自体が、そういう比較を推奨しないというメッセージにも思える。<br />
<br />
# もっとも、ポアソン分布の場合は今回のような選択肢があるけれど、二項分布の場合はGLMMを使うくらいしか対処の方法がない(怪しげなQAICの使用を回避するならば)。<br />
<br />
<span style="color: #3d85c6;"># 2014.06.24 追記:実際のデータのモデル化を考えた時、固定効果のみで推定値を求めるときGLMM (poisson)とglm.nbは異なる推定値になってしまうのは気になるが、もはやモデルの違いと割り切るしかない気がしている。双方ともGLMMなので、誤差はランダム効果でオーバーフィットしてしまうため、AICなどで比較して差異があっても実質的なモデルの良し悪しを反映しているのだろうか?現実のデータの場合は真の答えは分からないし、</span><span style="color: #3d85c6;">モデルを使用している以上は仮定を捨てられないので、</span><span style="color: #3d85c6;">悩んでも仕方がないと振り返って思います。</span><br />
<br />
# 参照:<br />
# lmerはresidual devianceに基づいたAICを返しているので、直接比較をするにはdevianceに基づいた計算にし直す必要がある<br />
# deviance: -2*logLik()<br />
# residual deviance: deviance()<br />
# AIC(glm())の場合は、deviance + 2*推定するパラメータ数<br />
# AIC(lmer())では、residual deviance + 2*推定するパラメータ数<br />
<br />
<br />
<br />
<br />
<br />
# 上記の数値実験を実行するためのプログラムはこちら:<br />
# 平均 exp(alpha + beta*x)とし、真の値はalpha = 2、beta = 0.5、あと負の二項分布の歪みのパラメータθ= 1/4、正規乱数ノイズの分散4も設定<br />
# (ただし、θや分散の値については今回はチェックはしない)<br />
<br />
library(MASS) # glm.nbの呼び出し<br />
library(lme4) # lmerの呼び出し<br />
<br />
OD <- function(alpha, beta, siml, distribution) {<br />
# 切片、傾き、シミュレーション回数、データ分布(1:ポアソン分布、2:負の二項分布、3:ポアソン分布+正規ノイズ)<br />
estim <- aic.d <- aic.rd <- prob.d <- prob.rd <- numeric(0) # 後でデータをくっつけるためのイントロンみたいなもの<br />
set.seed(36) # 乱数の種を指定(他の値にすると推定エラーが出るかもですglm.nbやlmerの問題)<br />
for (s in 1:siml) {<br />
Noise <- rnorm(100, 0, 4^0.5) # ノイズの分散を4とした<br />
x <- rep(c(0:9), each=10)<br />
y <- rbind(rpois(100, lambda=exp(alpha + beta*x)), # ポアソン分布<br />
rnegbin(100, mu=exp(alpha + beta*x), theta=1/4), # 負の二項分布<br />
rpois(100, lambda=exp(alpha + beta*x + Noise))) # 普通のポアソン分布にノイズを追加<br />
y <- y[distribution,] # データ分布の選択<br />
ID <- 1:length(y) # データの数だけ個体差(ID)を設定<br />
D <- data.frame(x, y, ID)<br />
pois <- glm(y ~ x, family=poisson, data=D) # GLM (poisson)<br />
negb <- glm.nb(y ~ x, data=D) # GLM (負の二項分布)<br />
pois.mm <- lmer(y ~ x + (1|ID), family=poisson, data=D) # GLMM<br />
estim2 <- c(cbind(coef(pois), coef(negb), fixef(pois.mm)))<br />
estim <- rbind(estim, estim2) # simlの数だけ推定値を重ねていく<br />
aic.p <- AIC(pois) + 2*logLik(pois) + deviance(pois) # residual devianceに基づくGLMのAIC<br />
aic.n <- AIC(negb) + 2*logLik(negb) + deviance(negb) # residual devianceに基づくGLM.nbのAIC<br />
aic.m <- -2*logLik(pois)[1] - deviance(pois) + AIC(pois.mm) # devianceに基づくGLMMのAIC<br />
Aic <- c(AIC(pois), AIC(negb), aic.m) # deviance に基づくAIC<br />
Aic2 <- c(aic.p, aic.n, AIC(pois.mm)) # residual deviance に基づく AIC<br />
aic.d <- rbind(aic.d, Aic)<br />
aic.rd <- rbind(aic.rd, Aic2)<br />
prob.d <- c(prob.d, which.min(Aic))<br />
prob.rd <- c(prob.rd, which.min(Aic2))<br />
}<br />
aic.d <- apply(aic.d, 2, mean)<br />
aic.rd <- apply(aic.rd, 2, mean)<br />
prob.d <- factor(prob.d, levels=factor(c(1:3))) # モデルを示す番号をカテゴリー化<br />
prob.d <- table(prob.d)/siml # AICモデル選択の結果を集計<br />
prob.rd <- factor(prob.rd, levels=factor(c(1:3))) # モデルを示す番号をカテゴリー化<br />
prob.rd <- table(prob.rd)/siml # AICモデル選択の結果を集計<br />
<br />
select <- rbind(aic.d, prob.d, aic.rd, prob.rd)<br />
colnames(select) <- c("glm.pois", "glm.nb", "glmm.pois")<br />
par(mfcol=c(1,2))<br />
plot(density(estim[,1]), lwd=4, xlim=c(-2,6), ylim=c(0,1), main="alpha (intercept)")<br />
for (i in 2:3) { lines(density(estim[,(2*i-1)]), lwd=4, col=i)<br />
abline(v=alpha) } # alpha値の図示<br />
plot(density(estim[,2]), lwd=4, xlim=c(0,1), ylim=c(0,6), main="beta (coefficient)")<br />
for (j in 2:3) { lines(density(estim[,(2*j)]), lwd=4, col=j)<br />
abline(v=beta) } # beta値の図示<br />
return(select) }Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-38318065187386165852014-02-10T14:03:00.003+09:002014-02-27T16:29:25.220+09:00負の二項分布のパラメータ推定:R と WinBUGSの比較過分散した整数データのモデリングによく使う負の二項分布、Rではglm.nb関数でできるが、WinBUGSの例をほとんど見ない。ちょうど使う用事があったので使い方をテストしてみた。<br />
<br />
# サンプルデータは、平均 exp(2)、dispersion parameter 1.5 の負の二項分布から抽出した N = 1000 のデータ(後で対数リンク関数を使用するのでexp()で表記)<br />
<span style="color: lime;">require(MASS)</span><br />
<span style="color: lime;">Y <- rnegbin(1000, mu=exp(2), theta=1.5) </span><br />
# 注:乱数なので実行する度に結果は少しずつ変わる<br />
<div>
<br />
<br />
<span style="color: #3d85c6;"># まずはRのglm.nb関数による推定</span><br />
<span style="color: lime;">require(MASS)</span><br />
<span style="color: lime;">summary(glm.nb(Y ~ 1))</span><br />
<div>
<br /></div>
<div>
<div>
Coefficients:</div>
<div>
Estimate Std. Error z value Pr(>|z|) </div>
<div>
(Intercept) 1.97519 0.02771 71.28 <2e-16 ***</div>
<div>
---</div>
<div>
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1</div>
<div>
<br /></div>
<div>
(Dispersion parameter for Negative Binomial(1.5895) family taken to be 1)</div>
</div>
<div>
<br /></div>
<div>
<span style="color: #3d85c6;"># 平均 1.97519(expの中身に相当)、dispersion parameter 1.5895 と順当な推定</span></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<span style="color: #3d85c6;"># 次にWinBUGS、負の二項分布だからdnegbin</span>、と試してみるが…(thetaがdispersion parameterに相当しそう)</div>
<div>
<div>
<span style="color: lime;">data <- list(N=1000, Y=Y)</span></div>
<div>
<span style="color: lime;">parm <- list(theta=rlnorm(1), alpha=rnorm(1))</span></div>
<div>
<span style="color: lime;"><br /></span></div>
<div>
<span style="color: lime;">model <- function() {</span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>theta ~ dgamma(1E-1, 1E-2)</span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>alpha ~ dnorm(0, 1E-6)</span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>for (i in 1:N) {</span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>Y[i] ~ dnegbin(p[i], theta) </span># 負の二項分布では平均をそのまま使えない</div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>p[i] <- theta/(theta+mu[i]) </span># pと平均muの関係(pは0~1の値を取る)</div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>log(mu[i]) <- alpha </span># muは正の値</div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>}</span></div>
<div>
<span style="color: lime;">}</span></div>
<div>
<span style="color: lime;">mcmc <- wbugs(data, parm, model, n.iter=10000, n.burnin=5000, debug=T)</span> # <a href="http://nhkuma.blogspot.jp/p/r2winbugs-macwinewindows-wbugs.html">自作ラッパー関数wbugsを使用</a></div>
<div>
<br /></div>
<div>
# order of negative binomial *** must be an integer というエラー</div>
<div>
# theta は負の二項分布の狭義の定義通り、整数でないといけないようだ</div>
</div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<span style="color: #3d85c6;"># 戸惑いつつも多項分布の事前分布で thetaが正の整数しか取らないようにし再計算</span></div>
<div>
<div>
<span style="color: lime;">data <- list(N=1000, Y=Y, q=rep(1,20)/20)</span></div>
<div>
<span style="color: lime;">parm <- list(theta=sample(c(1:5), 1), alpha=rnorm(1))</span></div>
<div>
<span style="color: lime;"><br /></span></div>
<div>
<span style="color: lime;">model <- function() {</span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>theta ~ dcat(q[]) </span># 事前分布1~20の整数(多項分布使用)</div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>alpha ~ dnorm(0, 1E-6)</span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>for (i in 1:N) {</span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>Y[i] ~ dnegbin(p[i], theta)</span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>p[i] <- theta/(theta+mu[i]) </span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>log(mu[i]) <- alpha </span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>}</span></div>
<div>
<span style="color: lime;">}</span></div>
<div>
<span style="color: lime;">mcmc <- wbugs(data, parm, model, n.iter=10000, n.burnin=5000)</span></div>
</div>
<div>
<br /></div>
<div>
<div>
mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff</div>
<div>
theta 2.000 0.000 2.000 2.000 2.000 2.000 2.000 1 1</div>
<div>
alpha 1.976 0.026 1.927 1.958 1.976 1.993 2.026 1 1500</div>
<div>
deviance 6032.409 1.446 6031.000 6032.000 6032.000 6033.000 6036.000 1 1500</div>
</div>
<div>
<br /></div>
<div>
<span style="color: #3d85c6;"># 平均 alphaの推定は1.976と正しく推定できているのだが、thetaは当然のように整数値となった…</span></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
# Rの?rnbinomなどを読み返してみると、"dispersion parameter (the shape parameter of the gamma mixing distribution). Must be strictly positive, need not be integer" とある。</div>
<div>
# gamma mixing distributionのshapeパラメータに相当するということか。</div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<span style="color: #3d85c6;"># なので、ポアソン分布にガンマ分布を混ぜた階層ベイズとして負の二項分布の推定を試みる</span></div>
<div>
<div>
<div>
<span style="color: lime;">data <- list(N=1000, Y=Y)</span></div>
<div>
</div>
<span style="color: lime;">parm <- list(shape=rlnorm(1,log(1),log(2)), alpha=rnorm(1), lambda=rlnorm(1000,log(1),log(2)),F) </span></div>
<div>
<span style="color: lime;">model <- function() {</span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>shape ~ dgamma(1E-1, 1E-2)</span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>alpha ~ dnorm(0, 1E-6)</span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>for (i in 1:N) {</span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>Y[i] ~ dpois(lambda[i]) </span># 負の二項分布を階層ベイズで表現</div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>lambda[i] ~ dgamma(shape, rate[i])</span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>rate[i] <- shape / mu[i]</span> # rate と shape、平均 mu の関係を表す</div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>log(mu[i]) <- alpha </span></div>
<div>
<span style="color: lime;"><span class="Apple-tab-span" style="white-space: pre;"> </span>}</span></div>
<div>
<span style="color: lime;">}</span></div>
<div>
<span style="color: lime;">mcmc <- wbugs(data, parm, model, n.iter=10000, n.burnin=5000)</span></div>
</div>
<div>
<br /></div>
<div>
<div>
mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff</div>
<div>
shape 1.587 0.087 1.422 1.522 1.587 1.647 1.755 1.002 1200</div>
<div>
alpha 1.975 0.028 1.917 1.957 1.975 1.993 2.029 1.003 1200</div>
</div>
<div>
# lambdaは省略</div>
<div>
<br /></div>
<div>
<span style="color: #3d85c6;"># 平均 1.975、shape 1.587 と 元の値やglm.nbの推定とほぼ一致、これが正解のようだ</span></div>
<div>
<br /></div>
<div>
# それにしてもややこしい…負の二項分布には異なる定義が多すぎる!</div>
<div>
<br /></div>
</div>
Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0tag:blogger.com,1999:blog-6412094354437530509.post-46089337167133868232013-09-30T15:36:00.001+09:002013-09-30T17:39:02.838+09:00混合分布のパラメータ推定:その4(複数種類のデータで分ける)<br />
(このページは URL を2013ベントス・プランクトン学会大会の要旨に掲載したものからリンクしています)<br />
<br />
# 4部構成にしてあります、ここではグループ数の推定がうまくいかない場合の対処について。<br />
# より基本形の例は、その1の方でまとめてあります。<br />
<a href="http://nhkuma.blogspot.jp/2013/08/blog-post.html"># その1:正規混合モデル(オーソドックスな基本形)</a><br />
<a href="http://nhkuma.blogspot.jp/2013/08/blog-post_13.html"># その2:ガンマ混合モデル(体長の増加に伴い分散も増加するモデル)</a><br />
<a href="http://nhkuma.blogspot.jp/2013/09/blog-post.html"># その3:前後関係からユーザがグループ数を仮定して正規混合モデルを実行</a><br />
# その4:体長と体幅など複数種類のデータで分ける(有効な最終手段)<br />
<br />
*******************************************************<br />
<br />
その1〜3では、ひとつの基準(体長とか、個体サイズを表す情報)を元にした分け方を紹介した。しかし、ここまでの方法だけではどうしても分けられないケースは出てくるだろう。<br />
<br />
ところで、<span style="color: #6fa8dc;">複数の基準…体長と体幅とか、体長と頭幅とか、体長と体重、のように1サンプルにつき複数の情報が得られるならば、より高い精度でコホート分離が可能である。こういうデータをお持ちの方 or 取ることが可能な方はぜひ試したらよいと思う。</span><br />
<br />
そもそも、ここまでで紹介した mclust は本来はそういう複数の基準で混合分布を解くのが基本のパッケージである。基準同士の相関関係も利用して分離してくれるので強力である。<br />
なお、多変量の混合分布なので、多変量正規混合分布(multivariate normal mixture distribution)のモデリングとなる。<br />
<br />
<br />
<br />
### まずはサンプルデータの設定<br />
<span style="color: lime;">N <- 100</span><br />
<span style="color: lime;">mean1.2 <- c(10, 1.0)</span><br />
<span style="color: lime;">mean2.2 <- c(13, 1.3)</span><br />
<span style="color: lime;">var1.2 <- matrix(c(1.2, 1, 1, 1.2), 2, 2)^2 </span><br />
<span style="color: lime;">var2.2 <- matrix(c(1.4, 1, 1, 1.4), 2, 2)^2</span><br />
<br />
### 母集団からのサンプリング、N1, N2 ずつ取ってきて混ぜた<br />
<span style="color: lime;">ratio <- c(0.8, 0.2)</span> # 8 : 2 の比率でサンプリングした<br />
# つまりトータル100で80、20個体のコホートが混じっている<br />
<span style="color: lime;">N1 <- N*ratio[1]</span><br />
<span style="color: lime;"> N2 <- N*ratio[2]</span><br />
<span style="color: lime;"><br /></span>
<span style="color: lime;">set.seed(1) </span># 乱数のタネを指定<br />
<span style="color: lime;">library(MASS)</span><br />
<span style="color: lime;">Length2 <- rbind(mvrnorm(N1, mean1.2, var1.2), mvrnorm(N2, mean2.2, var2.2))</span> # 1, 2を混ぜたサンプルを作る<br />
<br />
# (体長と体幅のようなサイズ関係)<br />
<span style="color: lime;">par(mfrow=c(1,2))</span><br />
<span style="color: lime;">hist(Length2[,1], col="gray")</span><br />
<span style="color: lime;">hist(Length2[,2], col="gray")</span><br />
<div>
<br /></div>
<div>
# 2峰であることは想像がつくが、ひどい混ざり方である。このケースではどちらか一方のみを用いた解析では分離できなかった。</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhVX66no4IwV_91CZynFi41cssIMU-GUdWX5baUviJ5dajUW7N41We9JFAv1YWp2eBeQuNnxbOKOtBZjCK7JKagQZ0e90bUzr3V7fg4p-D0s9WKEMcZwozfkOrRMykifuN5fVsERmfJ7GY/s1600/Rplot.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="181" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhVX66no4IwV_91CZynFi41cssIMU-GUdWX5baUviJ5dajUW7N41We9JFAv1YWp2eBeQuNnxbOKOtBZjCK7JKagQZ0e90bUzr3V7fg4p-D0s9WKEMcZwozfkOrRMykifuN5fVsERmfJ7GY/s320/Rplot.png" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
# 基本的なやり方はここまでと同じだが、</div>
<div class="separator" style="clear: both; text-align: left;">
# データは2基準のデータがセットであり2列になるようにする(上記のLength2のように)</div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
# BICのモデルはより複雑になる</div>
<div class="separator" style="clear: both; text-align: left;">
</div>
<div class="separator" style="clear: both;">
<span style="color: lime;">library(mclust) </span></div>
<div class="separator" style="clear: both;">
<span style="color: lime;">plot(mclustBIC(Length2)) </span></div>
<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgS9khmfZ-SWkhz7Lzc1RCKIjAMvHCE6PhyphenhyphentCLIzOLpE72w9LX_RVlULIHIOMgeua6gKAHIGMS_Wu0-XOmU7twOEZDfcoQWXsoBSaz6uIVsl3U25E0_Be5oc4ZI5oF_pq7mDWOI4DH5Nug/s1600/bic.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgS9khmfZ-SWkhz7Lzc1RCKIjAMvHCE6PhyphenhyphentCLIzOLpE72w9LX_RVlULIHIOMgeua6gKAHIGMS_Wu0-XOmU7twOEZDfcoQWXsoBSaz6uIVsl3U25E0_Be5oc4ZI5oF_pq7mDWOI4DH5Nug/s320/bic.png" width="304" /></a></div>
<div>
BICが最大となっているモデルのグループ分け数は2となり、ちゃんと分けられた。</div>
<div>
EEVなどのモデルはコホートの二次元分布のサイズ、形状、方向の仮定の組み合わせで表されている(詳しくは、<span style="color: lime;">citation("mclust")</span> で出てくるFraley et al. 2012のp.8のテーブルを参照)。</div>
<div>
<br /></div>
<div>
<br /></div>
<div>
# 推定値の見方もこれまでと同じ</div>
<div>
<div>
<span style="color: lime;">Mc <- densityMclust(Length2)</span></div>
<div>
<span style="color: lime;">mc <- summary(Mc, parameters=T)</span> </div>
<div>
<span style="color: lime;">mc</span></div>
<div>
<br /></div>
<div>
<div>
-------------------------------------------------------</div>
<div>
Density estimation via Gaussian finite mixture modeling </div>
<div>
-------------------------------------------------------</div>
<div>
<br /></div>
<div>
Mclust EEV (ellipsoidal, equal volume and shape) model with 2 components:</div>
<div>
<br /></div>
<div>
log.likelihood n df BIC ICL</div>
<div>
-309.8967 100 9 -661.2399 -664.1</div>
<div>
<br /></div>
<div>
Clustering table:</div>
<div>
1 2 </div>
<div>
86 14 # 正解は80と20なので少しグループ1の方が多目の判定</div>
<div>
<br /></div>
<div>
Mixing probabilities:</div>
<div>
1 2 </div>
<div>
0.8566939 0.1433061 </div>
<div>
<br /></div>
<div>
Means:</div>
<div>
[,1] [,2]</div>
<div>
[1,] 10.304959 14.246681 # 値の大きい方の基準は10、13が正解</div>
<div>
[2,] 1.126679 1.721643 # 小さい方の基準は1.0、1.3が正解なので、ちょっと粗い</div>
<div>
<br /></div>
<div>
Variances:</div>
<div>
[,,1]</div>
<div>
[,1] [,2]</div>
<div>
[1,] 1.7080561 0.9243498 # 1.7..と0.96...は2つの基準のそれぞれ分散</div>
<div>
[2,] 0.9243498 0.9646183 # 0.92..の2つは共分散</div>
<div>
[,,2]</div>
<div>
[,1] [,2]</div>
<div>
[1,] 1.1210113 0.9727447</div>
<div>
[2,] 0.9727447 1.5516630</div>
</div>
<div>
<br /></div>
<div>
# 推定は多少粗いものの、このデータでこれだけ分けられれば十分だろう。</div>
<div>
<br /></div>
<div>
# 図示の仕方も同じ、出てくるものは標高図になる</div>
<div>
<span style="color: lime;">plot(Mc, data=Length2)</span></div>
</div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjuaYy1MNJwjnqNY3HpbNqZZKY-kXR18PT_y7enEHKolX3WVtr1R1gXT3cYFtt-_rnZE0LWy3T84E0Y5tlJu6fnqO5UlhEf__coGn8zhyphenhyphen8Pq4mnUXcmsoM3dTwhElQISxsD0tGsaAaf3Tw/s1600/Rplot.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjuaYy1MNJwjnqNY3HpbNqZZKY-kXR18PT_y7enEHKolX3WVtr1R1gXT3cYFtt-_rnZE0LWy3T84E0Y5tlJu6fnqO5UlhEf__coGn8zhyphenhyphen8Pq4mnUXcmsoM3dTwhElQISxsD0tGsaAaf3Tw/s320/Rplot.png" width="310" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
横軸、縦軸に2つの基準を取り、二次元のヒストグラムを標高図にしている</div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div>
# ベイズ版は省略。複数基準ならMclustで十分だと思う</div>
<div>
<br /></div>
<br />Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com4tag:blogger.com,1999:blog-6412094354437530509.post-44344148685203932712013-09-25T16:29:00.000+09:002013-09-30T15:47:42.080+09:00混合分布のパラメータ推定:その3(グループ数をユーザ側で指定する場合のコード…)(このページは URL を2013ベントス・プランクトン学会大会の要旨に掲載したものからリンクしています)<br />
<br />
# 4部構成にしてあります、ここではグループ数の推定がうまくいかない場合の対処について。<br />
# より基本形の例は、その1、2の方でまとめてあります。<br />
<a href="http://nhkuma.blogspot.jp/2013/08/blog-post.html"># その1:正規混合モデル(オーソドックスな基本形)</a><br />
<a href="http://nhkuma.blogspot.jp/2013/08/blog-post_13.html"># その2:ガンマ混合モデル(体長の増加に伴い分散も増加するモデル)</a><br />
# その3:前後関係からユーザがグループ数を仮定して正規混合モデルを実行<br />
<a href="http://nhkuma.blogspot.jp/2013/09/blog-post_30.html"># その4:体長と体幅など複数種類のデータで分ける(有効な最終手段)</a><br />
<br />
*******************************************************<br />
<br />
(その1からの続きです)<br />
体サイズ頻度分布のヒストグラムでよくあることですが、高齢のコホートが他と比べてかなり少ない場合には、解析すると隣のグループに含められてしまうことが多々ある(想定するグループ数と異なると感じる場合…それは主観的かもしれないが)。<br />
<br />
<br />
もし、前の時点でのグループ数などの事前情報によって、そこに確かにコホートが存在するはずという仮定ができるのならば、グループ数は既知であるとして解析してもよいだろう。完全に仮定無しでの解析はなかなか難しい…混合分布問題は完全解決には至っていないと思う。<br />
<br />
<br />
その2でやり方を書いたような、グループ数指定(G=◯の項を加える)をしても、途中にもう一つグループを作られて、肝心の最高齢集団が無視されたりとなかなか言うことを聞いてくれない場合がある。<br />
<br />
その2のガンマ混合モデルならば解決できるだろうが、ベイズを使うのは敷居が高いだろうし計算時間もかかる。ここでは mclust(正規混合モデル)を用いた別の解決方法を紹介する。<br />
<br />
mclust (正規混合モデル)を使用して、グループ数をユーザ側で決めた場合の混合分布モデリングを行う。一番若い齢(体長などの値の小さい)のコホートから順に推定し、推定し終わったら除いて残りのデータで解析するというステップ・バイ・ステップ法で処理することによって、集団サイズの小さい最高齢のコホートも検出することができる。<br />
<br />
<div style="border: 2px dotted #CCCCCC; padding: 10px;">
### 母集団に用いるサンプルデータ<br />
<br />
<span style="color: lime;">N <- 200</span> # サンプル全体のN数<br />
<span style="color: lime;">mean1 <- 10</span><br />
<span style="color: lime;">mean2 <- 15</span><br />
<span style="color: lime;">mean3 <- 20</span><br />
<span style="color: lime;">sd1 <- 1.0</span><br />
<span style="color: lime;">sd2 <- 1.2</span><br />
<span style="color: lime;">sd3 <- 1.4</span><br />
<br />
### 母集団からのサンプリング、N1, N2, N3 ずつ取ってきて混ぜた<br />
<span style="color: lime;">ratio <- c(0.5, 0.3, 0.2)</span> # 5:3:2 の比率でサンプリングした<br />
<span style="color: lime;">N1 <- N*ratio[1]</span><br />
<span style="color: lime;"> N2 <- N*ratio[2]</span><br />
<span style="color: lime;"> N3 <- N*ratio[3]</span><br />
<span style="color: lime;"><br /></span>
<span style="color: lime;">set.seed(1) </span># 乱数のタネを指定<br />
<span style="color: lime;">Length3 <- c(rnorm(N1, mean1, sd1), rnorm(N2, mean2, sd2), rnorm(N3, mean3, sd3))</span> # 1, 2, 3を混ぜたサンプル</div>
<br />
<br />
<br />
<span style="color: #6fa8dc;">### mclust による正規混合モデリングを最若齢〜最高齢のコホートへと順に1つずつ推定するプログラム</span>(いったん読み込ませたら、Rを終了するまで何度でも再利用できる)<br />
<div style="border: 2px dotted #CCCCCC; padding: 10px;">
<span style="color: lime;">step.mclust <- function(data, n.group, var="V", priors=NULL) {</span><br />
<span style="color: lime;">require(mclust)</span><br />
<span style="color: lime;">x.min <- 5*floor(0.2*min(data))</span><br />
<span style="color: lime;">x.max <- 5*ceiling(0.2*max(data))</span><br />
<span style="color: lime;">x <- seq(x.min, x.max, by=(x.max - x.min)/100)</span><br />
<span style="color: lime;">summ <- numeric(0)</span><br />
<span style="color: lime;">for (g in 1:n.group-1) {</span><br />
<span style="color: lime;">mc <- summary(Mclust(data, G=n.group - g, modelNames=var), parameters=T)</span><br />
<span style="color: lime;">n <- round(mc$pro[1]*length(data)) </span># グループ1の推定個体数<br />
<span style="color: lime;">n2 <- round(mc$pro[2]*length(data)) </span># グループ2の推定個体数<br />
<span style="color: lime;">m <- mc$mean[1] </span># グループ1の推定平均値<br />
<span style="color: lime;">m2 <- mc$mean[2]</span><br />
<span style="color: lime;">s <- sqrt(mc$variance[1])</span> # グループ1の推定標準偏差<br />
<span style="color: lime;">s2 <- sqrt(mc$variance[2])</span><br />
# 1の中では確率密度が小さく、かつ2の方が確率密度が大きいものは2である<br />
# 1の中では確率密度が小さいが、2よりも確率密度が大きいものも1である<br />
<span style="color: lime;">d <- dnorm(data, m, s)</span> # dataがグループ1に属する確率密度<br />
<span style="color: lime;">d2 <- dnorm(data, m2, s2)</span> # 同、グループ2に…<br />
<span style="color: lime;">dd <- d > d2</span> # グループ2よりもグループ1に属する確率が高い個体<br />
<span style="color: lime;">p.rank <- rank(1/d)</span> # 確率密度の高い順<br />
<span style="color: lime;">gr <- data[dd ==T & p.rank <= n]</span> # グループ1に属する可能性が高い個体(でもまだグループ1が残ってる)<br />
# 面倒でも二段階抽出をしないとグループ1の左の裾が落ちてしまう<br />
<span style="color: lime;">gr2 <- setdiff(data, gr)</span> # dataからgrを除いた残り(この中にグループ1が残っている)<br />
<span style="color: lime;">r.remain <- rank(1/dnorm(gr2, m, s))</span> # 残り物の中で、グループ1に属する確率が高いものの順<br />
<span style="color: lime;">remain <- gr2[r.remain <= n - length(gr)]</span> # グループ1の残党<br />
<span style="color: lime;">gr <- c(gr, remain)</span> # remainをグループ1に追加<br />
<span style="color: lime;">data <- setdiff(gr2, remain)</span> # グループ2からremainを除外したものを次のループのdataにする<br />
<span style="color: lime;">summ <- rbind(summ, data.frame(n, m, s)) }</span><br />
<span style="color: lime;">rownames(summ) <- 1:n.group</span><br />
<span style="color: lime;">return(summ) }</span></div>
<br />
<div>
<br />
<div style="border: 2px dotted #CCCCCC; padding: 10px;">
### 使用法<br />
<span style="color: lime;">step.mclust(Length3, 3)</span> # データ、グループ数を指定するだけ<br />
<br />
### 出力結果(n:各グループの推定 N数、m:同 平均値、s:同 標準偏差)<br />
<br />
# n m s<br />
# 1 99 10.09472 0.8822265<br />
# 2 60 14.68063 1.4391196<br />
# 3 41 19.91919 2.3426576</div>
<div>
<br />
# 最後のグループの標準偏差の推定が良くない。<br />
# また、もっと最後のグループの集団サイズを減らすとけっこう推定は悪くなる…。とくに集団サイズの推定がよくない。<br />
そういうケースにもきっちりと対応したい場合には、その2のガンマ混合モデルでがんばるなど試行錯誤するしかないだろう。<br />
<br />
<br />
<br /></div>
</div>
Naoki H. KUMAGAI(熊谷直喜)http://www.blogger.com/profile/02435765594079260417noreply@blogger.com0