State Space Model

Summary 状態空間モデル(State Space Modeling)とは、モデル構造として状態過程と観察過程をもつ統計モデルの総称。JAGSとRスクリプトの備忘録。 State Space Modeling 生態学の時系列データ解析では、状態空間モデルがとても役立つ。なぜなら、「状態過程(State process)」と「観察過程(Observation process)」という二つプロセスを明示的に切り分けたモデル構造を持つからである。状態過程は、状態の時間変化をモデル化する(状態とは、集団サイズの増減のような実際には直接観察できないものを指す)。一方、観察過程では、状態に何らかの誤差が加わり、その結果としてデータが得られるという過程をモデル化する。例えば、状態過程では集団動態を表すモデルを当てはめ、そこに何らかの観察誤差が加わった結果として手元のデータが得られているという仮定を置く。これにより、集団増加率やその年変動など、将来予測に大きな影響を及ぼすパラメーターをより適切に推定できる(バイアスが少ない)。 JAGS script 次のような15年分の時系列データ(魚の個体数)があるとする。これに対し、状態空間モデルを当てはめてみる。 plot(dat$X, type = "o", ylab = "Fish abundance", xlab = "Year") State process:状態過程では、集団動態を指数増加関数としてモデル化することにする(対象に応じて変える)。各年の集団増加率\(r_t\)は、正規分布に従って年変動すると仮定し、その平均は\(\mu_r\)、分散は\(\sigma^2_r\)とする(対数スケールであることに注意)。 \[ ln~n_{t+1} = r_{t} + ln~n_{t}\\ r_t \sim Normal(\mu_{r},\sigma^2_r) \] Observation process:この集団動態モデルから予測される値に対し、観察誤差が加わった結果として観察値\(X_t\)が得られるとする。ここでは、観察誤差は過分散を考慮したポアソン分布(Overdispersed-Poisson distribution)に従うと仮定する。 \[ X_t \sim Poisson(n_{obs,t})\\ ln~n_{obs,t} \sim Normal(ln~n_t, \sigma^2_{obs})\\ \] JAGSのスクリプトで書き下す。 model{ # Prior mu.r ~ dnorm(0, … Continue reading State Space Model

恥ずかしい英語の間違い

Summary アメリカに来て二年が経ったが、いまだに恥ずかしい英語の間違いを繰り返してばかり。そんな失敗談のまとめ。 No thanks だれも知り合いがいない国際学会に一人で参加し、Burger Kingに行った時の話。当時(博士課程2年)、英語がおぼつかないながらも何とかハンバーガーを注文し、席で普通に食べていた。そうすると、店員の人がテーブルを回りながらお客に声をかけている。 "How XXXing?" 断片的にしか聞き取れていないのでなんのやり取りなのかさっぱりわからない。そしてついに私のテーブルにも来て、何やら同じようなことをしゃべりだす。何を言っているかわからない。だけど何か返さなければならない。そんな中、私は必死になにか手掛かりとなる情報を探し、その店員さんがケチャップとマスタードを持っていることに気が付き、これは「ケチャップいる?」とかその辺を聞いているに違いない!と合点した。味はしっかり濃い目だったので、調味料の類はいらないと思い、自身満々でこう答えた。 "No thanks" 周りの人たちが「ぎょっ」とし、店内の空気が凍ったのを覚えている。店員はバツの悪い笑顔を浮かべて去り、なにかこちらを睨みつけながら他の店員と話している。いったい、なにが起こったのだ。 あとあと思い返すと、おそらくあれは "How is everything?" と話しており、「ハンバーガーどう?」みたいに聞いていたのだと思う。それに対して私は、「全然ないわ」みたいな返しをしたことになる。ああ。 What are you working on? ポスドクの飲み会があり、それに参加して話していた時の話。目の前でビールを飲むブラジル人のポスドクがいたので、なんの研究をしているのか聞いてみようと思い、"What are you working on?"と聞いた。そうすると、 "I'm working on beer." と返ってきた。確かにそうだ。彼はビールを飲んでいる。無意味に現在進行形にしてはならないことを学んだ。 Everything 同じ部屋にいるポスドクにペンを借りようと思い、声をかけた時の話。貸して、というと(パンパンのペン立てを見せながら)どのペンがいいと言ってきたので、どれでもいいよ(Anything)、というつもりで "Everything" と答えてしまった。今でも恥ずかしい。 XXXgraXXX アメリカに住んで1年半がたち、次のポジション探しをしていた時の話。海外学振の任期は2年間なので、そろそろ次の当てを見つけなければならない。そう思い、アメリカと日本のポジションに応募していた。アメリカのポジションはすべてOnline Applicationなので楽なのだが、日本のポジションは海外からでも郵送で応募しなければならない。その時、応募書類のデータをCDに入れて送ってください、とするポジションがあり、(それだったらメールで出させてくれよ、とか思いながら)USPSの事務所に郵送手続きに向かった。 受付の姉ちゃんに封筒の中身は何だと聞かれCDだと答えると、「~でないこと」を証明する項目に署名してくれといわれるが、~の部分がXXXgrXXX contentsとしか聞こえず、「え?え?」と聞き返しまくった。そうすると、姉ちゃんの機嫌はどんどん悪くなってゆく。そこで目を落として書類を見ると Pornographic と、書いてある。あぁ、そうか。エッチなDVDを国外輸送されると困るから、そうではないことを証明してくれといっていたのか。どうやら私は、USPSの受付の姉ちゃんにPornographicを連呼させるという離れ業をやってのけたらしい。ごめんなさい。 Can I borrow a restroom? 友人の家に招かれ鍋をしていた時の話。トイレに行きたくなったので「トイレ貸して」のつもりで "Can I borrow a restroom?" というと、友人は軽く笑いながら"No"といった。え、マジ、と困惑していると、理由を教えてくれた。Borrowはペンなどのように「持ち手から”取って”借りる」場合を指すらしく、今回の場合のように使うと「トイレ丸ごと持ち出して借りてもいい?」という意味になるらしい。みなさん、Useを使いましょう。

Varying Intercept and Slope Model

Summary JAGSによるVarying Intercept and Slope Modelの書き方の備忘録。 Varying Intercept and Slope Model 切片と偏回帰係数がグループごとに変化するモデルを指す。書き下すと、次のような数式になる。 \[ y_i \sim Normal(\mu_i, \sigma^2)\\ \mu_i = \beta_{j(i),1} + \beta_{j(i),2}x_i\\ \beta_{k,j} \sim MVN(\mu_{\beta_k}, \sum) \] 切片および偏回帰係数である\(\beta_k\)がグループ単位\(j\)で変化すると仮定している。これにより、グループごとに切片および偏回帰係数の値が異なる(グループごとに応答変数\(y_i\)の平均水準と説明変数\(x_i\)の効果が異なる)ことを表現する。切片と偏回帰係数には多変量正規分布を事前分布とすることで、切片と偏回帰係数の相関も考慮する。ここでは次のテストデータに、このVarying Intercept and Slope Modelを当てはめてみる。 head(dat, 15) ## Y X groupID ## 1 1.6533548 3.715236 1 ## 2 1.9452796 1.334316 1 ## 3 1.0739305 7.079826 1 ## 4 … Continue reading Varying Intercept and Slope Model

Varying Intercept Model

Summary JAGSによるVarying Intercept Modelの書き方の備忘録。 Varying Intercept Model 切片がグループごとに変化するモデルを指す。書き下すと、次のような数式になる。 \[ y_i \sim Normal(\mu_i, \sigma^2)\\ \mu_i = \beta_1 + \beta_2x_i + R_{j(i)}\\ R_j \sim Normal(0, \sigma_R^2) \] 切片である\(\beta_1\)に対してグループ単位\(j\)で変化する変数\(R_j\)が加わった構造になっている。これにより、グループごとに切片の値が異なる(グループごとに応答変数\(y_i\)の平均水準が異なる)ことを表現する。意味としては同じだが、次のように書くこともできる。 \[ y_i \sim Normal(\mu_i, \sigma^2)\\ \mu_i = \beta_{1,j(i)} + \beta_2x_i\\ \beta_{1,j} \sim Normal(\mu_{\beta_1}, \sigma_{\beta_1}^2) \] ここで、切片のグループ間のばらつきを表す変数を、正規分布に従う確率変数とする。こうすることで、グループごとにそれぞれ切片を推定するのではなく、全体切片(\(\beta_1\)もしくは\(\mu_{\beta_1}\))と分散(\(\sigma_R^2\)もしくは\(\sigma_{\beta_1}^2\))という二つのパラメーターの推定で済ませることができる(グループ数が多くなるほどパラメーター数を節約できる)。ここでは次のテストデータに、このVarying Intercept Model(この場合GLMMと同じ)を当てはめてみる。 head(dat, 15) ## Y X groupID ## 1 0.7620270 5.1462547 1 … Continue reading Varying Intercept Model

J1 Visa: Two Year Ruleの落とし穴

Summary 注意事項 以下の記載内容は、筆者の個人的な経験をまとめたものにすぎません。間違いがないよう気をつけましたが、本記事の情報をもとに不利益を被ったとしても一切の責任を負いかねますのでご留意願います。 経緯 現在J1ビザでアメリカに滞在している私は、8月から始まる新しい職に就くためにH-1Bビザに切り替える必要があった。しかし、この関係で非常にややこしい問題にぶち当たった。J1ビザに伴う帰国義務(2年ルール)に課されていないものと信じていたのだが、実は課されていることが判明し、H-1Bビザに申請できない状態に陥ってしまったのである。この経緯について、今後同じような問題に陥る人が少しでも減るよう情報をまとめておきたいと思う。 J1ビザの帰国義務(2年ルール) J1ビザとは、知識の国際交流を目的に創設されたアメリカのビザのカテゴリーであり、その取得の容易さから多くの研究者がアメリカに長期滞在する際に利用している。しかし、その容易さと引き換えに制約がある。それが今回の問題の発端となった2年ルール(Two-Year Home-Country Physical Presence Requirement)である(参照)。このルールは、アメリカで経験を積んだのちに最低2年間は母国で過ごさない限り(もしくは帰国義務免除が承認されない限り)、Hビザやグリーンカードなど、特定のビザに申請できなくなるというものだ。ただし、このルールはJ1ビザで渡航している人すべてに課されているわけではなく、以下の条件のいずれかを満たした人に課されるとされている(アメリカ国務省のWEBページにより;参照)。 Government funded Exchange Program - You participated in a program funded in whole or in part by a U.S. government agency, your home country’s government, or an international organization that received funding from the U.S. government or your home country’s government. Specialized Knowledge or … Continue reading J1 Visa: Two Year Ruleの落とし穴

Regression Example in JAGS

Summary 様々な確率分布を仮定した単回帰モデルのJAGSスクリプトの例。基本的には、リンク先のスクリプトのLiklihoodの箇所を変更し、各種パラメーターの事前分布を指定してやればよい。 Gaussian Yのとりうる値の範囲:実数 \[ y_i \sim Normal(\mu_i, \sigma^2)\\ \mu_i = \beta_1 + \beta_2x_i \] Y[i] ~ dnorm(mu[i], tau) # tau is precision mu[i] <- b[1] + b[2]*X[i] Poisson Yのとりうる値の範囲:0以上の整数 \[ y_i \sim Poisson(\lambda_i)\\ ln\lambda_i = \beta_1 + \beta_2x_i \] Y[i] ~ dpois(lambda[i]) log(lambda[i]) <- b[1] + b[2]*X[i] Binomial Yのとりうる値の範囲:0以上N以下の整数 \[ y_i \sim Binomial(N_i, p_i)\\ … Continue reading Regression Example in JAGS

How to use JAGS: Simple Linear Regression

Summary JAGSの基本的な使い方の備忘録。ここでは、もっとも単純な単回帰モデルをJAGSで実装する。実装の手順は大きく二つに分けられる。(1)JAGSモデルの記述:JAGS言語を使ってJAGSのモデルファイルを記述する。(2)Rによる実装:RからJAGSのモデルファイルを走らせる。 JAGSモデルの記述 サンプルデータ ここではもっとも単純な例として、(GLMでも分析できる)単回帰モデルの分析をJAGSでする。ここで使用するサンプルデータは次のようになっている。 head(dat1) ## Y X ## 1 11.859762 6.635249 ## 2 4.682198 2.290394 ## 3 14.809464 7.288998 ## 4 7.549894 3.507948 ## 5 8.841815 3.150085 ## 6 13.473219 8.793556 plot(Y ~ X, dat1) このデータに対し、応答変数\(y_i\)を説明変数\(x_i\)で説明する簡単な統計モデルを考える(添え字\(i\)はサンプルのID番号)。応答変数\(y_i\)の誤差は正規分布に従うと仮定し、以下の統計モデルを考える。 \[ y_i \sim Normal(\mu_i, \sigma^2)\\ \mu_i = \beta_1 + \beta_2 x_i \] ここで、\(\beta_1\)は切片(intercept)を、\(\beta_2\)は回帰係数(slope)を表しており、これらが推定対象となる主なパラメーターである。 JAGSの記法 上の統計モデルをJAGS言語で書き下すと以下のようになる。 # JAGS … Continue reading How to use JAGS: Simple Linear Regression

Line

Summary 散布図にGLMなどの予測値を加える方法の備忘録。plot関数で図示したのち、abline関数やlines関数で上書きする。 Gaussian function “abline” First argumentに切片、Second argumentに傾きを入れる 正規分布を仮定したモデルの場合のみ有効 fit <- lm(Y ~ X, dat1) plot(Y ~ X, dat1) # estimated intercept and slope fit$coefficients ## (Intercept) X ## 0.07065496 0.20422560 # first argument: intercept, second argument: slope abline(fit$coefficients[1], fit$coefficients[2]) function “lines” 説明変数Xの最小値から最大値までの範囲の値を100等分(もしくはより多く)する変数xを作る 新しい変数xを基にしたdata.frameを作る predict関数で予測値をえる lines関数で書き加える fit <- lm(Y ~ X, dat1) # set new … Continue reading Line

Plot

Summary plot関数の使い方に関する備忘録 ポイント ポイントタイプを変える plot関数のpchを変更する plot(Y ~ X, dat1, pch = 1) plot(Y ~ X, dat1, pch = 19) ポイントサイズを変える plot関数のcexを変更する(デフォルトを1としたときの相対値を指定) plot(Y ~ X, dat1, cex = 3) フチの色を変える plot関数のcolを変更する plot(Y ~ X, dat1, pch = 21, col = "red") 中身の色を変える plot関数のbgを変更する(ただし、pch21など中身の色を変えられるシンボルを指定) #with border plot(Y ~ X, dat1, pch = 21, bg = "red") #with … Continue reading Plot

Ordinal Logistic Model

Summary 順序はあるが、定量的ではないデータに適用される統計モデル。例えば、環境状態を良い順に1-5のスコアで評価したデータなどが考えられる。詳細はこちらで解説されている。 function “polr” パッケージ:MASS ランダム効果:不可 その他:なし # sample data plot(as.numeric(dat1$Y) ~ X, dat1) # fitting with "polr" fit <- polr(Y ~ X, data = dat1) summary(fit) ## ## Re-fitting to get Hessian ## Call: ## polr(formula = Y ~ X, data = dat1) ## ## Coefficients: ## Value Std. Error t value ## X … Continue reading Ordinal Logistic Model