約 14,767 件
https://w.atwiki.jp/hayatoiijima/pages/20.html
Rパッケージ ここでは自分がこれまでにインストールしたパッケージを載せています。PCがクラッシュしたときのためです。 aod aodが含まれる。これも一般化線形混合モデルのときに使う。familyでnbinomとbetabin(成功,失敗)が使える。 car TypeIIやTypeIIIの分散分析が行える関数Anova()を含む。 coda lmer()の結果のオブジェクトを投入すると、経験ベイズ法による結果を算出してくれるmcmcsamp()を含む。 date 日付データを扱うのに必要な関数がそろっている。 e1071 skewness(歪度)やkurtosis(尖度)が計算できるようになる。実はscatterplot3dを含んでいる。 faraway Extending the linear model with Rの作者がオリジナルで定義した関数、データなど。逆ロジット変換(ilogit)が意外と便利だったりする。 geoRglm 空間情報を考慮したGLMMをおこなうのに必要。 glmmML glmmMLが含まれる。一般化線形混合モデルができるようになる。familyはbinomialとpoissonだけ。 gplots barplot2などが含まれる。エラーバー付き棒グラフの作成が簡単に。 Hmisc Rで得られた結果をLaTeX形式で出力するためのlatex()を含む。 labdsv 植生データ解析で使う。 lme4 lme4などが含まれる。lme4は階層を持ったエラーを持つ場合のGLMM(一般化線形混合モデル)を計算させることができる。 R2WINBUGS RからWinBUGSを動かすためのパッケージ。 scatterplot3d scatterplot3dが含まれる。3次元散布図が作れるようになる。実はe1071に含まれている。 survival 生存時間解析には欠かせない。 smatr II型の回帰分析、および関数関係の比較を行うパッケージ。 vegan これもlabdsv同様、植生データ解析で使う。 xtable Rで得られた結果をLaTeX形式で出力するためのxtable()を含む。
https://w.atwiki.jp/hayatoiijima/pages/41.html
このページは、こちらに移転しました。
https://w.atwiki.jp/hayatoiijima/pages/50.html
Rで作図 このページでは、Rでの作図方法について掲載します。データハンドリングともかぶるところがあるかも。 更新履歴:とりあえず移転 2011年2月13日からのアクセス数: - R でのグラフの作り方は、(一つの画面に何枚の図を描くかの指定)→ 高水準作図→ 低水準作図という順番に行います。 よく使う作図関数 hist() ヒストグラム。一変量の分布の形を見るのに使う。 d - read.csv("data.csv") hist(d$omosa) #omosaのヒストグラムが表示される。 boxplot() 箱ひげ図。カテゴリーごとの一変量の分布を見るのに使う。カテゴリーは何層にも階層化できる。 boxplot(ha ~ syori, d) #(描画する値 ~ カテゴリー, データフレーム名) #二重カテゴリーでもできる例を示す d$CO2 - rep(c("A", "B"), 112) boxplot(ha ~ syori + CO2, d) plot() 散布図。二変数の関係を見るのに使う。 plot(ha ~ ne, d) #カテゴリーごとに色やシンボルを変える plot(ha ~ ne, pch=as.numeric(syori), + col=as.numeric(syori), d) #pchはシンボルを数字で指定、colは数字や色名で色を指定。 #数字と指定できるシンボル、色の対応を見るためには、 library(Hmisc) show.pch() #warnings()が出ますが無視していいです show.col() #事前に色を指定しておくことも可能 #カテゴリーがアルファベット順に処理されることを利用。 iro - c("black", "red", "blue") plot(miki ~ omosa, xlim=c(0, 7), ylim=c(0, 2), + pch=as.numeric(d$syori)-1, + col=iro[as.numeric(d$syori)], d) #iro[*] とすると、iro の*番目に入れた色が使われます。 #as.numeric() で各カテゴリーを数字化して指定。 xyplot() カテゴリー別の散布図。library(lattice)に収録。 library(lattice) xyplot(ha ~ ne | syori, d) #(formula | カテゴリー名, データフレーム名) xyplot(ha ~ ne | syori + CO2, d) #当然カテゴリーは階層化できる。 pairs() データフレーム全ての列の組み合わせの散布図。 d - read.csv("data.csv") pairs(d) #この場合も色やシンボルをカテゴリーごとにすることができる。 pairs(d, pch=as.numeric(d$syori), + col=as.numeric(d$syori)) 作った図の保存の仕方 ゐんどーづの場合:作った図の上で右クリックして、ビットマップにコピーを選び、コピーする。そして、ペイントなどに貼り付けて、ファイルとして保存する。 コマンドラインから指示する場合:dev.copy()を使う。 PDFファイルとして出力する場合。 #作図する dev.copy(pdf, file="好きなファイルネーム.pdf") dev.off() #作図後すぐに行うのがコツ PNGファイル(GIFとほぼ同等。ただしより高機能) #作図する dev.copy(png, file="好きなファイルネーム.png", + width=500, height=500) #widthとheightは適当に。 dev.off() #作図後すぐに行うのがコツ ちなみに、作った図のフォントのタイプを変えたい場合、斜体か、太字か、ゴシックかぐらいは作図のときに指定できるが、もう少し細かい指定は、 #作図する dev.copy(pdf, file="好きなファイルネーム.pdf", + family="Times") dev.off() というように、family=""の中に指定をすると変更できる。ちなみに、これはPDFファイルの場合のみ有効な指定らしい。 また、PDFでは特に明示しない場合、縦横比が1 1の図が生成されるが、これを変えたい場合は縦と横の相対的な長さの比率をwidthとheightで与えてやる。 #作図する dev.copy(pdf, file="好きなファイルネーム.pdf", + family="Times", width=10, height=5) dev.off()
https://w.atwiki.jp/hayatoiijima/pages/49.html
Rと統計解析 このページでは、Rでの統計解析を中心に解説します。 更新履歴:ページを作成(2011/1/29) 2011年1月29日からのアクセス数: - Rと統計解析 GLMMこんなときはGLMMだ! RにおけるGLMMの実行関数lmer()の推定結果は要注意 各種一般検定関数こんなときは検定だ! t検定 U検定対応のあるt検定 ANOVA(F検定) Kruskal Wallis検定 多重比較 ANOVA(N 元配置の分散分析) 相関 回帰 非線形回帰 各種クロス集計表の検定 GLMM GLMについては知っているものとして解説します。 こんなときはGLMMだ! 影響を検討したい変数以外に、応答変数に明らかに影響を与えている要因がある。 その要因に関して、要因として評価はしたくないが、無視もできない。 このような、「興味はないが応答変数に影響を与える変数」のことを、Random effectといい、これを含んだGLMをGLMM(Generalized linear mixed model)と呼びます。私たちが取るほとんどのデータは、実はRandom effectを含んでいます(林学関係でこれがないと断言できる状況なんてないと思うが)。 Random effectの例~世の中Random effectだらけ! 個体差(同じ処理区なのに......) ブロック間差 その場所に固有な何か 時間(年)変動 ということで、GLMMはいまや解析の最も基本的な部分を占めつつあります。 RにおけるGLMMの実行関数 glmmML() パッケージglmmMLに収録。library(glmmML)で使える。 長所 stepAIC()など既存の関数が適用しやすい。安定している。 短所 扱える誤差構造はbinomial、poissonのみ。 入れられるRandom effectは1個だけ 基本形 #REはRandom effect、dはデータフレーム名 glmmML(y ~ x1 + ..., cluster=RE, family=**, d) lmer() パッケージlme4に収録。library(lme4)で使える。 長所 扱える誤差構造がgaussian、binomial、poissonなど多い。 Random effectは複数指定可能(しかも切片と傾きどちらも可)。 短所 関数の仕様が特殊なため、既存の関数が適用できないことが多い。 基本形 lmer(y ~ x1 + ... + (1 | RE), + family=**, method="Laplace", d) それぞれの解析例を示します。 library(glmmML) res - glmmML(surv ~ light, cluster=FL, + family=binomial, d) summary(res) coef se(coef) z Pr( |z|) (Intercept) -2.369 0.376 -6.300 2.99e-10 light 13.111 2.126 6.166 6.99e-10 Standard deviation in mixing distribution 1.2e-05 Std. Error 0.2354 Residual deviance 256.1 on 221 degrees of freedom AIC 262.1 GLMと大きな差はありません。Null devianceが出力されないぐらいが違いでしょうか。 library(lme4) res - lmer(surv ~ light + (1 | FL), + family=binomial, method="Laplace", d) res Generalized linear mixed model fit using Laplace Formula surv ~ light + (1 | FL) Data d Family binomial(logit link) AIC BIC logLik deviance 262.1 272.4 -128.1 256.1 Random effects Groups Name Variance Std.Dev. FL (Intercept) 5e-10 2.2361e-05 number of obs 224, groups FL, 29 Estimated scale (compare to 1 ) 0.9923228 Fixed effects Estimate Std. Error z value Pr( |z|) (Intercept) -2.3773 0.3767 -6.311 2.77e-10 *** light 13.1629 2.1299 6.180 6.41e-10 *** --- ...... 一部省略 こちらは大分様子が違います。まずモデルの当てはまりを示すAICなどが示されており、その下にRandom effectがどの程度のばらつきを説明しているかが表示され、その下に説明変数に関する情報が載せられています。結果のオブジェクトをfixef()に入れると説明変数が、ranef()に入れるとRandom effectの値が表示されます。 lmer()の推定結果は要注意 2010年10月13日追記。East_scrofaさんのところで色々と実験がなされています。結果を用いるときは、妥当性をよく確かめて下さい。 http //d.hatena.ne.jp/East_Scrofa/20101012 http //d.hatena.ne.jp/East_Scrofa/20101009 各種一般検定関数 こんなときは検定だ! すでに検討したい要因が絞られていて、その要因の影響の差を検討したい。 厳密な制御環境にあり、検討したい要因以外の影響を考慮する必要がない。 指導者(あるいはEditor)に検定しろって言われて逆らえない^^;)(このケースが一番多いか?) 検定の原理に関しては逐一説明しませんが、基本的には 差を見ようとする集団間で差がないと仮定する。 両集団からランダムサンプリングしてきて差をとった値(あるいはもう少しこねくり回した値。いわゆる検定統計量)が既存の理論分布(正規分布またはそれに類するものとすることが多い)に従うと仮定する。 実際の測定値から検定統計量を計算したときに、両集団に差がないと仮定したときの検定統計量の分布においてどこに位置するかを調べ、すごくはしっこ(いわゆる95 %点よりも端。つまり珍しいゾーン)だったら、「両集団に差がないという仮定が間違っていた」と結論する。 ということです。 まずは以下の感じでオブジェクトを作っておく。 d - read.csv("data.csv") d2 - split(d, d$syori) t検定 t.test()を使う。 d0 - cbind(d2[[1]], d2[[2]]) t.test(omosa ~ syori, data=d0, var.equal=TRUE) var.equalをTRUEにしておけば通常のT検定、FALSEの場合(こちらがデフォルト)は、いわゆるWelchの検定となる。 U検定 wilcox.test()を使う。U検定は数学的にはWilcoxonの順位和検定と呼ばれる検定と同様であり、Rの中では、Wilcoxonの順位和検定が行われる。 wilcox.test(omosa ~ syori, data=d0) 対応のあるt検定 t.test()の中で、paired = TRUEとすることで実行できる。 ANOVA(F検定) aov()を使う。ただし、ANOVAは各群の母分散が等しい(等分散性)ことを仮定するので、それを確かめるために、bartlett.test()によって、等分散性の検定を事前に行っている。 bartlett.test(omosa ~ syori, d) result - aov(omosa ~ syori, d) summary(result) #結果の完全な取出しにはsummary() Kruskal Wallis検定 kruskal.test()を使う。 result - kruskal.test(omosa ~ syori, d) 多重比較 パラメトリックな多重比較法として最も一般的に使われるTukey法の場合。 result - aov(omosa ~ syori, d) TukeyHSD(result) ノンパラ検定の場合によく使われるStee-Dwass 検定の場合。決まった関数は用意されていないので、群馬大学の青木先生が作ったスクリプトを利用するのがよい。リンクはこちら。 これを一度R上にまるごとコピペしてEnterしておくと、以下のような方法でSteel-Dwass検定が行えるようになる。多重比較したいデータのある列を最初に、カテゴリーを示す列を次に持ってくる。カテゴリー変数は数値化しておく必要があるので、as.numeric() で囲っている。 (スクリプトをR上にコピペ) Steel.Dwass(d$miki, as.numeric(d$syori)) ANOVA(N 元配置の分散分析) library(car)に含まれるAnovaとlm()を使う。 library(car) result - Anova(lm(miki ~ omosa * syori, d), type="II") result 式の中の*は単独項と交互作用項を含める記号。 なら交互作用項のみが投入される。 相関 cor.test()を使う。 cor.test(d$omosa, d$miki) 回帰 lm()を使う。 result - lm(miki ~ omosa, d) summary(result) 非線形回帰 nls()を使う。用意したデータセットでは非線形回帰するようなシチュエーションがないので適当に以下のようにデータを作る。 Height - c(0, 0.3 12.3) Age - c(0, 11.8, 26.6, 32, 35.2, 37.6, 40.2, + 43, 46, 49.4, 53.4, 59.2, 71.4, 91) d - data.frame(Height, Age) plot(Height ~ Age, d) んで、 result - nls(Height ~ a/(1+b*exp(-c*Age)), + start=c(a=10, b=100, c=0.5), d) というようする。start=c()で、初期値を決めている。うまく収束しなかった場合は、何かしらのエラーメッセージが出る。 各種クロス集計表の検定 chisq.test()やfisher.test()を使う。 d - matrix(c(100, 20, 70, 50), ncol=2, byrow=T) rownames(d) - c(``Engineer , ``Agriculture ) colnames(d) - c(``Male , ``Female ) d #データの表示 chisq.test(d) 2£2 の表なら、fisher の正確確率検定の方がよい。 fisher.test(d)
https://w.atwiki.jp/hikkosi-hyouban/pages/19.html
アーク引越センター口コミ・評判 【総合評価】★★★ 2008年テレビ東京系「TVチャンピオン2」の引越し屋さん選手権で、厳しい主婦審査員15名 から圧倒的支持を受け、優勝した経歴を持ちます。これ、私も放送見ましたがさすがと思い ましたね。出場は松本引越センター、アリさんマークの引越社、サカイ引越センター、アーク引越センター、ダック引越センター、ハート引越センター、ベスト引越センターの7社だったんですが最後まで残ったのはアリさんとアーク。 結局主婦審査員の審査は9票VS6票でアークの勝利でした。ちなみに第1回、2回は松本引越 センターが優勝しているとか。 【口コミ・評判一部紹介】 ◎作業員がとても感じがよく、テキパキと手際よく進めてくれた ◎予定外の不用品の引き取りも気持ち良く対応してくれ、助かった ××作業開始時間が6時間も遅れた ××元バイトです。社員はちゃんとした教育を受けていない、ヤクザみたいな教養の低い人 しかいません。勤務初日、未経験のバイトに対しても罵倒・罵声の嵐。全く何の説明も されていないのに、「キルティング持ってこい!」「ジャバラさげる!」等訳のわから ないことを指示し、それが何か聞くと「役立たずが!」と一蹴されました。 プロの引越し評判・口コミ比較 アーク引越センター編
https://w.atwiki.jp/hayatoiijima/pages/39.html
輪読会Bayesian methods for ecology このページでは、2008-2009年に北大農学部で行われたMichael A. McCarthyの「Bayesian methods for ecology」の輪読会の様子を紹介します。今度は失敗しないぞ、と。 で、この本ではWinBUGS単体で使うことを念頭において書かれているのですが、いまどきWinBUGSを単体で使うことは(特殊な事情がない限り)まずないでしょう。ということで、我々はR2WinBUGS経由でこのゼミは進めます(本文には全く記載はありません、要注意!)。アップされている資料に関しても、R2WinBUGSを使っての解析用のコードになっています。 R2WinBUGSの使い方に関しては別ページ(こことかこことか)で解説しています。 このゼミは聞くだけでも構いませんし、途中参加も構いません。幅広い方のご参加をお待ちしております。 累積訪問者数(from 06/11/2008): - 今日来てくれた人: - 輪読会Bayesian methods for ecology ゼミの日時・場所目次と担当者 第1章Introduction条件付確率 Bayesの定理 Bayes統計の4要素 Example1 Example2 第2章Critiques of statistical methods 第3章Analysing averages and frequencies 第4章How good are the models? 第5章Regression and correlation 第6章Analysis of variance 第7章Mark-recapture analysis 第8章Effects of marking frogs 第9章Population dynamics 第10章Subjective priors ゼミの日時・場所 場所:北海道大学北方生物圏FSC2F大会議室 目次と担当者 Introduction(高橋) Critiques of statistical methods(平川) Analysing averages and frequencies(江口・渡辺) How good are the models?(深谷) Regression and correlation(宮田・川森) Analysis of variance(上野・飯島) Mark-recapture analysis(南波) Effects of marking frogs(赤坂) Population dynamics(志田) Subjective prior(森) Conclusions 第1章Introduction 担当者:高橋Ch1.pdf 条件付確率 :Cが起こる確率 :Dが起こる確率 :Dが起きた元でCが起こる確率 :CとDが両方起こる確率 とすると、 となり、これは と書き換えることができる。 Bayesの定理 Bayesの定理は、条件付確率を拡張したものである。条件付確率の定義に従って考えると、 であるが、はCとDが両方起きる確率であるから、 のように、CとDを逆にしても当然成り立つ。 とは等しいから、 が成り立ち、 となる。 ここで、Cはある仮説が正しいときに起こる事象(Ha)とし、Dはデータが発生する事象と仮定すると、 となり、Bayesの定理が導かれる。 言葉を再定義すると、 :事後確率 :事前確率 :尤度関数。仮説Haが正しいとした時にDが得られる確率。つまり、パラメータが真のときにデータが得られる確率なので尤度関数。 :周辺確率。データが得られる確率、というよくわからない確率。単にPr(Ha | D)を確率変数(積分して1にする)ための補正項として考えてもよい。 となる。 Bayes統計の4要素 例えば 朝、家を出る前にニュースで天気予報を聞く(prior) 玄関を出て空模様を見る(data) ごにょごにょ考える(model) 雨が降りそうと結論(posterior)し、傘を持って出る。 Example1 Example2 第2章Critiques of statistical methods 担当者:平川 Ch2.pdf box22A.Rbox22Amodel.txt box22B.Rbox22Bmodel.txt 全てまとめて落とすにはこちらch2.zip 第3章Analysing averages and frequencies 担当者 前半(~P77):渡辺 3章前半.ppt 後半(P77~):江口 ch3後半.lzh 第4章How good are the models? 担当者:深谷 Chapter4.lzh 第5章Regression and correlation 担当者 前半(~p141):川森 後半(p141~):宮田Chapter5 p141-.zip 第6章Analysis of variance 担当者 前半(p.158-170):上野Ch6_1.zip 後半(p.170-194):飯島Ch6_2.zip 第7章Mark-recapture analysis 担当者:南波Ch7.lzh 第8章Effects of marking frogs 担当者:永美 Ch8.zip 第9章Population dynamics 担当者:志田 Ch9.lzh 第10章Subjective priors 担当者:森 #ref error :ご指定のファイルが見つかりません。ファイル名を確認して、再度指定してください。 (Chapter10.lzh)
https://w.atwiki.jp/hayatoiijima/pages/17.html
統計解析ログ ここでは統計解析に関して学んだことを書き散らしていきます。正確さについては保証できませんのであしからず。 作成途中が多すぎるな...少しずつ直していきます。 GLM GLMM(作成途中) R作図小技集(内容を改変中) Rパッケージ(えんどれす?) Ⅱ型の回帰分析
https://w.atwiki.jp/hayatoiijima/pages/40.html
R2WinBUGSの使い方 このページは、R上からWinBUGSを操作することを可能にするパッケージ、R2WinBUGSの使い方を解説するページです。 モデル式の部分を加筆中 累積訪問者数(from 2008/12/17): - 今日来てくれた人: - R2WinBUGSの使い方 はじめに 作業手順 モデル式を書き下ろす 今回のデータ構造とその解析目的 モデル式の書き方 BUGSコードの基本的なルール 決定論的指定 確率論的指定 実装例 Rにデータを読ませる 初期値・パラメータの設定 bugs()の仕様 結果の確認 要約値の出力 収束診断図 簡単な事後分布図 chainごとの事後分布図 はじめに 計算機の進歩により、Bayes統計(というかMCMC計算)が可能になり、その利用に関する関心が急速に高まっています。しかし、実行方法の難しさが多くのユーザーへの拡大を妨げているといえるでしょう。そこでMCMC(Gibbsサンプリングのみですが.......)計算を行うためのWinBUGSが開発されました。WinBUGSは単体で結果の表示、作図などを行ってくれます。しかし、それら以外の結果が取り出したいケースもありますし、何よりRで使えたら論文などへの二次利用も容易となるでしょう。 そこで開発されたのがR2WinBUGSです。名前は本当にそのまんまと言う感じです(2はtoをもじったものでしょう)。R2WinBUGSを使うと、データはRで準備すればよく、WinBUGSでのMCMC計算の過程などをRに格納することが可能となり、利便性が格段に向上します。 といっても、R2WinBUGSはR2WinBUGSでいろいろお作法というか指定しなければならないことがあります。ここでは初心者が陥りがちな点を記述するつもりです。 作業手順 モデル式を書き下ろし、.bugという拡張子をつけて保存する。できた.bugファイルはRのワーキングディレクトリに保存する Rを立ち上げ、必要なデータをリスト形式で用意する 初期値の設定(指定しないこともできますが、複雑なモデルでは指定しないとまず収束しません)、結果として見たいパラメータを指定する library(R2WinBUGS)またはlibrary(arm)でパッケージを読み込む bugs()関数を使って計算を行う 収束程度などを作図により確認する なお、armというパッケージは、R2WinBUGSや作図で使うcodaなど、複数のパッケージの集まりです。個別にパッケージを導入するのが面倒くさい方は(私もその一人^^;)、armを導入するといいと思います。 各部分の詳細について以下で述べますが、せっかくですから実データを使って解析をしてみましょう。使うデータはこちら。R.csv モデル式を書き下ろす ここがある意味最も難しいところです。自分が考えているモデルを表現する必要がありますが、なれないうちは中々書くことができません。こればかりは慣れですので、頑張って習得してください。 今回のデータ構造とその解析目的 今回使うデータの解析の目的は、 植物に窒素を与える区と与えない区で、葉の窒素濃度に差が出るか? という問題です。窒素を与える区と与えない区でそれぞれ同種の植物が数個体あり、各個体から5枚ほど葉がサンプリングされています。 実際にファイルを開いてみると、データは以下のようになっていることがわかります。 No:個体のID Leaf:各個体の中での葉のID Ntreat:窒素処理の有無 Treatment:別の処理(今回は使いませんがあえて入れてあります。理由は後述) Ncon:葉の窒素濃度 モデル式の書き方 BUGSコードの基本的なルール BUGSコードを記述するファイルは、以下のように作る必要があります。 model { ......(モデルの中身) } このように記述したファイルを、.bugという拡張子をつけて保存します。 で、問題となるのはモデルの中身の書き方ですね。それを以下で解説します。 決定論的指定 要するに代入する部分です。 - で指定すると決定論的指定になります。決定論的指定にしては誤差などを挟む余地はないので、主に自分が考えているモデルを記述する部分になります。 今回の仮説は「植物に窒素を与える区と与えない区で、葉の窒素濃度に差が出るか?」でした。なので単純に考えられるモデルを表現すると、以下のようになりそうです。 Ncon[i] - a + coef * Ntreat[i] #aは切片(窒素が与えられないときの葉の窒素濃度)、coefは係数、iは各行を示す aとcoefだけ決まれば葉の窒素濃度がわかりそうです。でも、こんなに現象は単純でしょうか?例えば以下のようにして図を描いてみると、 plot(Ncon ~ Ntreat, col=No, d) #Noは各個体 Ntreatだけではっきりと区別できているわけではないですね。また、色の違いは個体の違いですが、個体によってずいぶんと窒素濃度が違うように見えます。 そうすると、まずモデルの予測値とデータの間にはある程度の誤差が含まれ、かつ個体によっても葉の窒素濃度に違いがありそうです。それを記述しようとすると、 Ncon[i] - mu[i] + error[i] mu[i] - a + coef * Ntreat[i] + catID[No[i]] #mu[i]はモデルの予測値、catIDは個体によって異なるカテゴリカル変数 となりそうです。 しかし、 -は代入ですから、未知の切片(a)や窒素処理の有無の係数(coef)、誤差(error)、個体間差(catID)は表現することができません。そこで、以下の確率論的指定が必要になります。 確率論的指定 確率分布を使う箇所は、確率論的指定をする必要があります。 ~ で指定すると確率論的指定になります。具体的には、 データに対し、モデルから得られる予測値を確率分布を使って当てはめる箇所 推定パラメータなど、確率分布から値を生成する箇所 といった場所は確率論的指定をする必要があります。 さて、これを踏まえて先ほどのモデルを考えると、aやcoef、error、catIDは確率論的指定を使えば表現できそうです。ただし、確率論的指定という名称からもわかるように、「どの確率分布」を使うかを決める必要があります。 特に、事前情報が無くどのような確率分布を使うかわからないときは、できるだけ無情報な確率分布を当てはめることが多いです。今回は、「平均が0、非常に分散が大きい正規分布」を用いることにします。 Ncon[i] ~ dnorm(mu[i], 0.001) #確率論的指定 mu[i] - a + coef * Ntreat[i] + catID[ID[i]] #決定論的指定 a ~ dnorm(0.0. 0.001) #確率論的指定 coef ~ dnorm(0.0, 0.001) #確率論的指定 catID[1] ~ dnorm(0.0, 0.001) #確率論的指定 ...... catID[21] ~ dnorm(0.0, 0.001) #確率論的指定 さてさて、いきなりレベルが上がったように感じると思いますが....... dnorm(mu, 1/σ)とは、平均mu、分散σの正規分布に従う値を生成する関数です。二項分布であればdbin()、ポアソン分布であればdpois()などがあります(詳しくはマニュアルを参照して下さい) 正規分布の分散については、逆数で指定するというルールがあります。なので、分散が0.001とは小さいように思えるか知れませんが、実際には1/0.001=1000です。 これらがわかったところで、もう一度上記のモデルを検討しましょう。気になるのは個体差の係数catIDですね。このままだと個体ごとに係数を独立に推定することになり、パラメータを21個も消費することになります。 しかし実際には個体差は興味のある対象ではないので、21個もパラメータを消費して個別に推定したいわけではありません。そこで、「個体差の係数は全個体でみると平均値0、分散σの正規分布に従う」と考えます。こうすると、正規分布に従うという制約が加わりますので、パラメータ数の消費を抑えることができます。 for (id in 1 21) { catID[id] ~ dnorm(0.0, tau) } またまたわけのわからない形に見えますが....... for (記号 in 数) {繰り返す内容}は、{}で囲まれた部分について「記号」に「数」を最初から最後まで代入する関数です。 つまり、idに1が入り、2が入り、......、最後に21が入ります。 さて、この21個の係数達の関係(似ているか似ていないか)は、係数の分散が小さいか大きいかという問題と言えます(分散が大きければ この分散については特に事前の情報がないので、やはりある確率分布に従う値として推定させます。ここでどのような確率分布を使うかですが、使われる分布は正規分布との積になるので、共役事前分布を使うと計算が楽になります。 共役事前分布とは、別の確率分布との組み合わせによっては、積の前後で同じ確率分布になる物です。よく知られているものとしては、 for (id 1 21) { catID[id] ~ dnorm(0.0, tau) } tau ~ dgamma(0.001, 0.001) 実装例 完成品がこちら。bayes.bug Rにデータを読ませる R上ではデータをデータフレーム形式で用意していることがほとんどかと思いますが、bugs()に渡すデータは以下の点に気をつけて準備してください。 データはリストの形で渡す必要があります 使用しないデータが含まれていてはいけません 文字列はデータとして使えません。カテゴリーデータであってもas.numeric()などで数字にしておいて下さい 初期値・パラメータの設定 後述しますが、WinBUGSが途中でこける場合の原因の多くが、初期値設定にあります。stochastic node(モデル式中で ~ で表されている部分)には初期値をきちんと与えましょう。 また、繰り返し回数やパラメータ数が多い時に全てのパラメータを結果としてみようとすると、bugs()計算が終了してRに戻るときに、「メモリーがいっぱいです」などと言われて結果が見れないことがあります。初期値はstochastic node全てに必ず与えないといけませんが、見るパラメータは絞りましょう。 bugs()の仕様 基本的には以下のbugs()関数を使ってMCMC計算を行います。bugs()はR2WinBUGSに含まれる関数です。 res - bugs(d, inits, para, "***.bug", n.chains=3, n.iter=3000, n.thin=10, n.burnin=1000, bugs.directory="c /自分がインストールしたディレクトリ", working.directory=getwd(), debug=TRUE) それぞれの項目の意味ですが、 d:データが入っているリスト inits:初期値を指定したオブジェクト para:結果としてみるパラメータを指定したオブジェクト ”***.bug”:モデル式を指定したbug ファイルの名前 n.chains:同時に走らせるシミュレーション数 n.iter:MCMC サンプリングの繰り返し数 n.thin:MCMC サンプリングを、何個おきに行うか n.burnin:MCMC サンプリングの内、最初のいくつを使用しないか bugs.directory:WinBUGS をインストールしたフォルダ(R2WinBUGS のデフォルト設定はc /Program files/WinBUGS14/なので、ほとんどの人はいじる必要がないはず。) working.directory:*.bugファイルのありかを示す。普通はRの作業フォルダにしてあるはずなので、getwd()としておけば、Rの作業ディレクトリを参照してくれる。新しい?R2WinBUGSではこの指定がないと動かない。 debug:TRUE にすると、MCMC計算が終了、あるいはうまく回らなかったときにWinBUGS が立ち上がったままになり、原因を特定できる(ことがある) 各種回数を設定するところが最も悩ましいところかと思います。どのぐらいで収束するかはデータとモデルのよしあしから決定されます。n.chainsは複数本とするのが普通で、データが大きいなど特別な理由がない限り1にはしないで下さい。普通は3ぐらいだと思います。 繰り返し数を増やせば収束しやすくなりますが、その分計算時間が増大します。今回のように少ないデータであれば問題となりませんが、1000を超えるようなデータではよく考えないと、計算時間が数時間に及びます。 n.thinはサンプリングの自己相関を除くために設定しますが、最終的に得るデータの数を減らす(きちんとシミュレーションした上で)という意味あいもあります。 慣れない内は上手くまわせないことが多いですから、debug=TRUEとしておきましょう。WinBUGSにサンプリングの図が出てくれば計算はきちんと終了していますので、WinBUGSを閉じるとRで作業できるようになります。一方、失敗した場合も、debug=TRUEとしておけばWinBUGSは立ち上がったままになりますから、原因を特定するのに役立ちます(見てもわからない場合も多々ありますが^^;)。 結果の確認 ちゃんと回ってよかったよかった、というわけにはいきません。収束していなければ意味がありません。収束を確認する最もよい方法は、MCMC計算の過程や事後分布を作図によって図示することです。 要約値の出力 たとえばresというオブジェクトにbugs()の結果を付与している場合、 res と打つと要約された推定値などが出力します。表示される桁数を変えたい場合は、 print(res, digits.summary=*) と打ちます。*の数字が出力桁数になります。 これを図で見る場合、 plot(res) とします。これを行うと、各パラメータの80%信頼区間(何故80%かはわかりませんが)、収束判定に使うR-hat値(古谷 2008を見てください)などが表示されます。 収束診断図 names(res) とすると、bugs()の結果にどのような情報が含まれているのかわかります。このうち、 sims.matrixには「chainの平均のサンプリング値」 sims.arrayには「chainごとのサンプリング値」 が含まれています。サンプリング図を示せば、収束したかは判断できそうですね。 簡単な事後分布図 こちらはchainの平均の値が使われていますから、厳密な収束診断に用いるべきではありませんが、簡単に描画することができます。 plot(as.mcmc(res$sims.matrix)) とすると、パラメータごとのサンプリング過程と事後分布図が表示されます。 as.mcmc()というのは、データをmcmcオブジェクトという形式に変換するための関数です。mcmcオブジェクトをplot()すると、自動的にサンプリング過程や事後分布図が表示される、という仕組みです。 chainごとの事後分布図 こちらはchainごとのサンプリング値が取り出せるので、収束診断として最も適当な図です。ただし、図にするには少々工夫が要ります。 chainごとのデータは配列(array。行列matrixは二次元だが、それを多次元に拡張したもの)で保存されています。これをplot()で読める形に変換します。 res$sims.arrayを表示させると、3次元のデータになっていることがわかります。ではそれぞれの次元は何を表しているかというと、 res$sims.array[サンプリング数, chain数, パラメータ数] となっています。 先ほどのmcmcオブジェクトに変換するための関数as.mcmc()は、一度に1chainのデータしか変換することができません。そこで、指定した列に関数を適用できるlapply()を使い、各chainのデータをmcmcオブジェクトに変換します。 lapply()は、 lapply(使いたいリスト形式のデータ, 関数) という書き方が基本です。ただし、今回のように、as.mcmc()を全部に適用するのではなく、各chainごとに、というように、自分で動作を指定する場合は、 lapply(関数を適用する列番号, function(x) .....(動作の中身)) という書き方をします。 lapply()は結果をリスト形式で返しますが、そのままでは描画することができません。そこで、mcmc.list()という関数を使ってmcmcリストというこれまた変な形式のリストに変換して、一つのリストとします。これをplot()に放り込めば望むような図が出力されます。 以上を踏まえ、 plot(mcmc.list(lapply(1 res$n.chains, #1からres$n.chainsまで(要は指定したchain数) function(x) as.mcmc(res$sims.array[, x, ]) ) ) ) とすると、chainごとに自動で色分けされたサンプリング過程と事後分布図が表示されます。
https://w.atwiki.jp/hayatoiijima/pages/36.html
GLMM ここでは、一般化線形混合モデル(GLMM)について解説します。 GLMMGLMM概論 GLMMを行うRパッケージglmmML特徴 Random effectの指定方法 aod特徴 Random effectの指定方法 lme4特徴 Random effectの指定方法 GLMM概論 GLMMを行うRパッケージ glmmML aod lme4 glmmML glmmMLパッケージ内にあるglmmML()という関数を使う。 特徴 familyはbinomialとpoissonが使える Random effectは一個しか指定できない Random effectの指定方法 result - glmmML(formula, family=**, cluster=**, data=**) ように書き、cluster=**に、Random effectの入ったラベルを入力する。 aod aodパッケージ内にあるnegbin()という関数を使う。 特徴 familyはbinomialとnegative binomialが使える。 Random effectの指定方法 負の二項分布の場合 result - negbin(formula, ~ **, data=**) ように書き、~ **に、Random effectの入ったラベルを入力する。 lme4 lme4パッケージ内にはlmer()という関数があり、これを用いてGLMMを行う。 特徴 familyはgaussian、binomial、poissonが使える。 Random effectを複数指定できる Random effectは、切片に対して、およびある説明変数の傾きに対して設定できる めちゃくちゃしんどい計算させるとRごと落ちる場合がある Random effectの指定方法 ランダム切片の場合 result - lmer(y ~ x1 + x2 + ... + (1|**), family=**, data=**) のように書き、(1|**)の**に、Random effectの入ったラベルを入力する。 ランダム傾きの場合 result - lmer(y ~ x1 + x2 + ... + (0 + x*|**), family=**, data=**) のように書き、(0 + x*|**)のx*には説明変数、**にはRandom effectの入ったラベルを入力する。
https://w.atwiki.jp/hayatoiijima/pages/51.html
Rでのデータハンドリング このページでは、Rでのデータの整形・加工・プログラミング(少しだけ)方法について説明します。 更新履歴:とりあえず移転(2011/2/13) 2011年2月13日からのアクセス数: - Rでのデータハンドリング データの構造1個のデータの種類 データの塊の種類 個々のデータの変換データの生成 種類の変換 NA処理 データフレームの整形・加工データフレーム同士の結合 データフレームの横または縦展開 プログラミング繰り返し命令 文字列処理 データの構造 1個のデータの種類 ちょっとイメージしにくいかも知れませんが、1個のデータを取ってみても種類がいくつかあります。具体的には以下のようです。 数字:numeric。1、3.5、-10など 文字列:character。"Tekito"、"Darui"など。 要因:factor。見た目は文字列っぽいのだが、水準(順位)が付いている。 論理値:logical。TRUEまたはFALSE。 欠損:NA。これが入っていると動作しない関数があるので要注意。 非数:NaN。1/0とかやっちゃったりすると発生する。やっぱり関数の動作を妨げるので要注意。 データの塊の種類 と、一つ一つのデータに種類があるわけですが、それらの集まりに関しても種類があります。具体的には以下のようです。 ベクトル:c()。構造などがない一つながりのデータ。あらゆる種類のデータを含めるが、同一ベクトル内で異なった形式の値を混在させることはできない。むりくり異なった形式の値を混在させると、一番ランクの低い形式に合わされる。例えば、文字列が入っている場合は数字も文字列扱いになってしまう。 test - c(1, 2, "Tekito", TRUE) test [1] "1" "2" "Tekito" "TRUE" test[1] + test[2] Error in test[1] + test[2] 1 + 2 [1] 3 行列:matrix()。行列(縦と横)の形があるデータ。数字以外は含めない。 データフレーム:data.frame()。見た目は行列だが、数字以外も含める。実用的には最も多用する形式。行や列のラベルあるいは番号で行や列を取り出せる。行数の違う列が混在することはできない。 リスト:list()。構造や並び方が違うデータも無理矢理一つにできる。自分で作ることはあんまり無いが、統計や作図関数の結果を返したオブジェクトはリストになっていることが多い。 とりあえずこうしてみると、現在読み込んでいるdはデータフレームですが、PlotやSpなど一つ一つの列は「ベクトル」と見なすことができます。 個々のデータの変換 データの生成 コロン( ):等差数列を生成する r***:乱数を発生させる。正規分布はrnorm(データ数, 平均値, 分散)、二項分布はrbinom(データ数, 上限値, 生起確率)、ポアソン分布はrpois(データ数, ラムダ)。他にも色々。 1 10 [1] 1 2 3 4 5 6 7 8 9 10 5 (-4) #負の数ももちろんオッケー [1] 5 4 3 2 1 0 -1 -2 -3 -4 1 10*10 #かけ算すれば、1以外の等差数列もできる [1] 10 20 30 40 50 60 70 80 90 100 1 + 1 10*2 #工夫すればいろんな等差数列ができる [1] 3 5 7 9 11 13 15 17 19 21 rnorm(100, 100, 10) #平均値100, 分散10の乱数を100個作る rbinom(100, 1, 0.8) #生起確率が80%で、0または1を取る事象(成功と失敗、裏と表とか)を100回行う rpois(100, 1) #ラムダ(平均値)が1のポアソン分布に従うデータを100個発生させる。 種類の変換 as.character():文字列に変換。要因として認識されてしまった数字の文字列を数字として認識させるときにも使う。 #ダミーデータの生成 No - as.factor(1 100*2) y - rnorm(100, 10, 5) d - data.frame(No, y) summary(d) No y 2 1 Min. -1.182 4 1 1st Qu. 6.782 6 1 Median 9.436 8 1 Mean 9.604 10 1 3rd Qu. 12.468 12 1 Max. 23.609 (Other) 94 #Noは要因として認識されている #これを単に数字にしようとすると、要因の「順番」の数字になってしまう as.numeric(d$No) #このような数字を「見た目のまま」数字に変換するには、 as.numeric(as.character(d$No)) as.numeric():数字に変換。作図でよく登場。 as.factor():要因に変換 NA処理 is.na():NAかどうかをTRUEまたはFALSEで回答 ifelse(条件式, TRUEの場合, FALSEの場合):条件式を与え、TRUEの場合とFALSE場合に異なる動作を行わせる。 データフレームの整形・加工 データフレーム同士の結合 cbind():データフレーム同士を横方向で結合 rbind():データフレーム同士を縦方向で結合 merge():2つのデータフレームに共通なデータに合わせる形で結合 データフレームの横または縦展開 reshape():データフレームを楯または横方向に展開する。具体例を見て下さい。 プログラミング 繰り返し命令 for():同じような動作を繰り返す時に使える。 split():データフレームをカテゴリーごとに 文字列処理 paste(, , sep=""):2つ以上の要素を、sep=で指定した文字列で区切って結合する(sep=""とすれば間に何も入らない)。 sub: substring(文字列, 最初の位置, 終わりの位置):文字列を左から指定した位置で切り取る