チラシの裏の落書き日記

統計とか,研究とか,日常についての備忘録的なもの。

Stanでマルチレベル(マニュアルコピペ)

参考:
https://mc-stan.org/docs/2_18/stan-users-guide/multivariate-hierarchical-priors-section.html
1.13 Multivariate Priors for Hierarchical Modelsより
高速なコードにはなっていないが,これで基本的にはOK。
添字が慣れたものではないが,読めばわかるような感じになっている。
cross classifiedモデルへの拡張や,適当なlink functionを仮定して従属変数が連続量以外でもOKな状況に持っていけるはず。

#
# Stanでマルチレベル lkj分布を使ったりなんだり
#
library("rstan")

# データ
library(haven)
mlmdata <- read_dta("https://stats.idre.ucla.edu/stat/examples/imm/imm10.dta")
write.csv(mlmdata,file="imm.csv",row.names = F)
dat <- read.csv(file = "imm.csv")
head(dat)

# データの用意
N <- nrow(dat)
x <- dat[,c("ses","homework","sex","white")]
K <- ncol(x)
J <- length(unique(dat$schid))
u <- dat[,c("meanses","public" )]
u <- u[!duplicated(u[,1]),]
L <- ncol(u)
y <- dat$math

uniq_sch_id <- unique(dat$schid)
jj <- rep(NA, N)
for (j in 1:J){
  jj[dat$schid==uniq_sch_id[j]] <- j
}
dat_stan <- list(N = N,
            K = K,
            J = J,
            L = L,
            jj = jj,
            x = x,
            u = u,
            y = y)
# int<lower=0> N;              // num individuals
# int<lower=1> K;              // num ind predictors
# int<lower=1> J;              // num groups
# int<lower=1> L;              // num group predictors
# int<lower=1,upper=J> jj[N];  // group for individual
# matrix[N, K] x;               // individual predictors
# row_vector[L] u[J];          // group predictors; uはJ×Lの行列でデータを指定する。
# vector[N] y;    
help(package="rstan")

cat("data {
  int<lower=0> N;              // num individuals
  int<lower=1> K;              // num ind predictors
  int<lower=1> J;              // num groups
  int<lower=1> L;              // num group predictors
  int<lower=1,upper=J> jj[N];  // group for individual
  matrix[N, K] x;               // individual predictors
  row_vector[L] u[J];          // group predictors
  vector[N] y;                 // outcomes
}
parameters {
  corr_matrix[K] Omega;        // prior correlation
  vector<lower=0>[K] tau;      // prior scale
  matrix[L, K] gamma;           // group coeffs
  vector[K] beta[J];           // indiv coeffs by group
  real<lower=0> sigma;         // prediction error scale
}
model {
  tau ~ cauchy(0, 2.5);
  Omega ~ lkj_corr(2);
  to_vector(gamma) ~ normal(0, 5);
  {
    row_vector[K] u_gamma[J];
    for (j in 1:J)
      u_gamma[j] = u[j] * gamma;
    beta ~ multi_normal(u_gamma, quad_form_diag(Omega, tau));
  }

//  for (n in 1:N)
//    y[n] ~ normal(x[n] * beta[jj[n]], sigma);
// Speed UP!!!
 {
    vector[N] x_beta_jj;
    for (n in 1:N)
      x_beta_jj[n] = x[n] * beta[jj[n]];
    y ~ normal(x_beta_jj, sigma);
  }
}
",file="multilevel.stan")

fit1 <- stan(
  file = "multilevel.stan",  # Stan program
  data = dat_stan,    # named list of data
  chains = 2,             # number of Markov chains
  warmup = 1000,          # number of warmup iterations per chain
  iter = 2000,            # total number of iterations per chain
  cores = 2,              # number of cores (using 2 just for the vignette)
  refresh = 1000         # show progress every 'refresh' iterations
)

stan_trace(fit1, pars=c("gamma"))
stan_trace(fit1, pars=c("Omega"))
stan_trace(fit1, pars=c("tau"))
stan_ac(fit1, pars=c("gamma"))
stan_ac(fit1, pars=c("Omega"))
stan_ac(fit1, pars=c("tau"))
stan_hist(fit1, pars=c("gamma"))
stan_hist(fit1, pars=c("Omega"))
stan_hist(fit1, pars=c("tau"))

print(fit1, pars=c("gamma", "Omega", "tau"), probs=c(.1,.5,.9))
plot(fit1, pars=c("gamma", "Omega"))

以下のチュートリアルも非常に参考になる
https://web.stanford.edu/class/psych201s/psych201s/papers/SorensenEtAl.pdf

とある論文を読んで考えたこと

下のリンクの論文を読んで,ベイズ統計について,少し気になる記述が多かったので,気になった点や考えたことをメモしておく。
https://www.jstage.jst.go.jp/article/psychono/37/2/37_37.25/_pdf

ベイズ統計学については東京工業大学の渡辺澄夫先生のWebページが有用。
以下のページのすべてのリンクを熟読しておきたい。
http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/joho-gakushu2.html
また,「ベイズ推論:いつも何度でも尋ねられること」のページに書いてあることは一日一度読み上げたい。
http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/nandodemo.html

5.まとめ

(1) ベイズも最尤も正しい統計的推論ではないが、どちらでも推論を行うことはできる。

(2) 正しい統計モデルや正しい事前分布は存在しないが、それでも推論を行うことはできる。

(3) 推論の結果は間違っているが、どのくらい間違っているかを調べる方法はある。

(4) 統計モデリングとは、適切な統計モデルや事前分布を設計することで 間違いの度合いを小さくしようと試みること。真実が見つかるわけではない。

それはそれとして,本文をみていく。


心理学は心という 目に見えないものを扱う 学問であるから ,測定における 誤差に確率分布を仮定して真値に近づくという 論拠が必要だった。

のっけから日本語がよくわからない。

一方, ベイズ統計学の礎であるベイズの定理は古く18世紀に明らかになったものであるが計算機の圧倒的な発展と 確率的プログラミング言語の拡充によって,様々な学問分野においての応用が期待されるようになったのは, この数年の話である 。

WinBUGSが出てきたのは,1996年あたり。(Spiegelhalter D., Thomas A., Best N., and Gilks W. 1996b. BUGS 0.5: Bayesian inference Using Gibbs Sampling–Manual (version ii). Medical Research Council Biostatistics Unit, Cambridge.)この数年というには古い気がする。拡充といっているから,他のソフトも想定している?

これまでの利用方法は頻度主義的統計学とも 呼ばれ,我々にとっての新しい統計学ベイズ統計学はそもそも考え方を異にするとされる 。我々はこれをどう受け止め, どのように付き 合っていけばよいのか。

統計学(or 統計的な解析方法?)のことを指しているのだろうが,より具体的にはなんの利用方法か明確でない。
また,我々とは誰だろうか。

本論文では, まず 2章でこれまでの心理統計の問題点を指摘し, これに対してベイジアンのアプローチが有利な特徴を持っていることを示す。

心理統計ー>心理統計学
また,心理統計学ベイジアンアプローチという比較しようがないものを比較しているのが気になる。

2.1. 従来の統計法の短所
心理学において再現性の問題が取りざたされる 以前から ,帰無仮説検定の利用の仕方に問題点があることは指摘されて来た。

なぜか,「従来の統計法」というものが,「帰無仮説検定」のみを指すかのように記述されている。というか,「従来」ってなに?

しかし事前になんの想定もしないことがすなわち , データ生成過程(尤度) をも 想定しないことにつながっており , このことが結果に大きな差異を生む。

データ生成過程と尤度を同一視するのはどうなのだろうか。また,データ生成過程をも想定しないでどのように統計的推測を行うのだろうか。
そのような想定をしていない統計家(or 分析者)が存在するのだろうか。

このように,同じデータからでも 結論が異なることに注意が必要である 。この三つの状況は,それぞれ想定する確率空間が異なるため,同じ 事象であっても , それが生じる 確率が異なる 。 ベイジアンは確率分布が異なるのだから 結果が異なるのは当然と 考えることができるが,帰無仮説検定の枠組みではデータ生成過程を取り 込めない。

分析の想定が違うので当たり前では?
この例では,以下のような研究の目的があり,それを検証するためにコイン投げを行っていることを想定している。

コイントスをして,表が出るか裏が出るかを記録し,このコインがいずれの結果も 公平で偏らないものであるかどうかを検証する例を考えてみよう 。

だとすれば,適切な想定は二項分布を想定することだろう。
Negative binomial distributionを想定したり,5分間コインを投げ続けるのは,この目的に対して何を明らかにしようとしているのかわからない。筆者は,「ベイジアンは確率分布が異なるのだから 結果が異なるのは当然と 考えることができるが,帰無仮説検定の枠組みではデータ生成過程を取り込めない。」と述べているが,何を対立させて議論しているのかわからない。ベイジアンでなくとも,結果が異なるのは当たり前に理解できる。

これまでは実験計画の中にその過程を統制する仕掛けを埋め込み,方法や手続きを詳述することで実験者の考えを推察するしかなかった。 ここでデータの生成過程である 確率分布を明示すれば,実験者がどのような背景でデータを得ているかをも 明示することができる 。

実験計画法は分散分析を行うためのデータ収集方法で,誤差の分布に正規分布を仮定している。ここに実験者の想定を推し量る営みなどないと感じられる。このことは後述の変数変換の部分とも関わってくる。さらに言えば,実験計画は系統誤差を偶然誤差に変換する巧妙な手続きによって,実験の効果を検証できるようにした鮮やかな手続きと考えられる。
また,想定している確率分布を明示したからといって,データを得る背景についてはなんの情報も得られない。データがどのように得られたのかは,実験手続きやmethodの詳細な記述から理解するのが普通の論文の読み方だと思われる。

しかしそれを明示する 必要がない状況では,方法で「このような前提のもとで行われた」 という 主張があればそれを受け入れるしかない。

これにも疑問が残る。どういった意図にせよ,特定の現象を検討するために不適切な方法が取られていたのであれば,批判や議論がきちんとなされるべきことであろう。

研究者の悪意を疑うわけではないが,頻度主義的なやり方では,意図的でない場合でも誤った結果に陥る可能性が高いと 言わざるを得ない。 この点が帰無仮説検定の最大の欠点である,と筆者は考えている 。

ベイジアンのほうがよいことへの伏線だと思うが,ベイジアンだからといって問題がないわけではない。また,誤った結果というのはなんだろうか?タイプⅠ・Ⅱエラーについては,検定論の枠組みでしっかりと議論がなされているし,方法論も存在する。ネイマンピアソンの検定の枠組みに則れば,タイプⅠエラーをコントロールした行動ができる。

ベイジアンの本質を短い言葉で表現すると , わからないことを確率で積極的に表現すること ,といえるだろう 。群間の差,効果の大きさなど,実験者にとって知りたい実質的な 違いを確率分布として 表現するのであって,仮想的な確率空間の議論をするのではないため, その解釈はより 単純である 。確率分布で表現するということは,確率分布のパラメータで表現することであり , どういった分布がどのようなパラメータをもつのか, という 基礎的な知識は必要であるが,一度それを把握してしまえば, あとは実質的な値の検証になる 。

標語的すぎてあまりよくわからない。実質的な違いを確率分布で表現するとは?ベイズ統計学を用いたとしても,分析に使う確率分布は”仮想的”(分析者が仮定するもの)であって,何かの実態が存在するわけではない。その意味で,ベイズ統計だろうが検定に用いられる帰無分布だろうが,人間が想定したものに変わりはない。なぜ片方は理解が単純で,片方は理解が難しいといえるのだろうか。また,「確率分布で表現するということは,確率分布のパラメータで表現することであり」とはどういうこと?順序統計量を用いたノンパラメトリック的な考え方もあるし,計算機的にブートストラップ的な考え方も存在する。さらにいえば,ノンパラメトリックベイズのように,パラメタが無限に存在している分析方法もある。こういった状況を想定はしていないのだろうか?また,「あとは実質的な値の検証になる 。」と述べているがどういうことだろうか。仮説検定であっても,確率分布を想定して,精密な分布が導出できていればそれを用いることができるし,場合によっては漸近論に依拠した検定を構成することもできる。

帰無仮説検定の仮定と 同じように, データが正規分布に従うとしよう 。

細かい話だが,データが他の確率分布に従っていたとしても,母平均の差を検討するための標本分布は構成できると思う。

これらを数式で表すなら ~

~以下は数式なのだろうか(私の数式感が狭い可能性はある)。また,暗にデータが独立同一の分布に従うとしているが,明示されていない。

ここでわからないものはμA, μB, σの三つであるから , これらに事前分布を置いてデータと組み合わせると , μA, μB, σの事後分布が得られる 。 この事後分布は, データを踏まえて考えた各群の平均値・分散の取りうる 値の確率を表しており ,実質的な大きさそのものである 。 もし …なら 確率pで云々, という 回りくどい表現ではないため,初学者にもわかりやすいし誤解が生まれにくい。

そもそもデータ生成過程自体だってわかっていないのに,わからないことはこれだけだ,というのはあまりよくわからない。また,「 μA, μB, σの事後分布」は,「 μA, μB, σの同時事後分布」としたほうがよいと思う。気になるのは,「実質的な大きさそのもの」といっているが,実質的とはなんだろうか。そして,「もし …なら 確率pで云々, という 回りくどい表現ではないため,」といっているが,この分布の仮定自体が「もし」なのであるため,ベイジアンのほうがわかりやすいと言われると違和感がある。条件付き確率という点では,仮説検定の帰無分布もベイズ統計の事後分布も同じだと思うのだが。

ベイジアンの考え方は, データを得ることで事前情報をアップデートしていくという 営みであるから , 研究の累積や再現性の問題にも 貢献する 。

研究の累積については,メタ分析があるし,再現性の問題は科学の構造的な問題のようにも感じられるので,ベイジアンだからいいというわけではないと思う。

第一に,例えば変数変換など, データを検定モデルに合わせる 必要がない。対数変換や角変換など,検定分布に合わせるためのテクニックがあったが,一般化線形モデルが登場したことで,最近は目にすることは少なくなっている 。 しかしそうした習わしは, そもそも 分布を自由に選ぶベイジアンにとっては不自然な振る 舞いでしかない。 データ 生成メカニズムを自然にデザインすればよいのである 。

変数変換については,結果の直接的な解釈が難しくなるという意味ではわかるが,データを上手く変換すれば,厳密によく知られた帰無分布を使うことができるというのは無視できない利点だとも考えられる。ここで仮定している一般化線形モデルも十分にデータ生成モデルを考えているモデルと考えることができる。2値データに対して,ロジスティック回帰分析をするか,プロビット回帰をするか,はたまた他の分布を仮定するか,ベイジアンでなくとも考えることであると思う。

第二に,検定論ででてくる 仮定や補正から 無縁でいられることが挙げられる 。 t検定の手続きは, データの正規性の検定を行い,分散の等質性の検定を行い, Welchの補正をかけるという 段階を踏まねばならない。 ベイジアンはどの事前分布と 尤度, データの組み合わせで事後
分布が得られるかについて,一つの定理を適用するだけでよい。 ある 場合は 5% より 大きければ喜び,続いてある 場合は 5% より 小さければ喜ぶといった仮定の組み合わせや, 補正という 手続きに悩むことがない。例えば先の例で,分散が等しいかどうかわからないのであれば,
Ai~N(μA, σA2), Bi~N(μB, σB2) とすればよいのである 。 σA,σBの事後分布を見て, どの程度異なるのか,等質とみなしてよいかどうかは, その研究領域の共通見解で判断すればよい。 あるいは, モデル適合度の観点で,異なる 分散を推定したモデルと 分散を同じとしたモデルの, どちらがデータに適合していたかを持って考えのヒントとすることもできるだろう 。

t検定についていえば,段階的に検定をおこなうこと自体がすでに誤りで,Welchの検定をはじめから行う必要がある。分散が等質でない2群の平均値差の推測はベーレンスフィッシャー問題といわれ,たしかにこんな簡単な設定でも困ったことが起こるというのは,問題だと思われる。ただし,「 σA, σBの事後分布を見て, どの程度異なるのか,等質とみなしてよいかどうかは, その研究領域の共通見解で判断すればよい」といっているが,唐突に数値的な分析を放棄したのはどうしてだろうか。ベイジアンは「ある 場合は 5% より 大きければ喜び,続いてある 場合は 5% より 小さければ喜ぶといった仮定の組み合わせや, 補正という 手続きに悩むことがない。」といっているが,事後分布の95%信用区間を使った意思決定を行っている人も多いように感じられるがどうなのだろうか。信頼区間が0を含むかどうかは同じ有意水準の検定と同等であるが,ベイジアンの事後信用区間についても同様と考えられるが。また,「あるいは, モデル適合度の観点で,異なる 分散を推定したモデルと 分散を同じとしたモデルの, どちらがデータに適合していたかを持って考えのヒントとすることもできるだろう 。」というのは,頻度論であっても十分に実行されてきた慣行だろう。情報量基準や,残差分析,尤度比検定などいくらでも方法はある。なぜ,これらを考えずに,あたかもベイズ統計学だけがデータ適合を考えてきたかのように記述を行っているのだろう。

第三の利点として,下位検定を事前に設計しておく 必要が生じないことが挙げられる 。

突っ込むのが疲れてきたが,実験計画の話とそうでない場合の検定のはなしとあらゆることがごっちゃになって議論がなされていると思う。実験データを解析するのに,ベイジアンな方法を使って,下位検定を想定しないで,好きな組み合わせだけを調べることが良いことなのかは私にはわからない。興味がある条件の組み合わせの検討方法として,対比分析なども存在するが,そういったことは考慮していないだろうか。

これに対してベイジアンでは,検証する 対象である 事後分布は常に一つである 。得られた事後分布,事後予測分布の, どの箇所をどのように切り 分けて(周辺化して) 比較しても , そのことで結論が変わることにはならない。

事後分布は一つではない。データ生成メカニズムや事後分布の想定によって変化する。また,周辺化して比較しても結論が変わらないと行っているが,結論とはなんだろうか。全体的に言葉がゆるふわであまり理解ができない。

もしどこかに交互作用が出るのではないかというような,探索的な仮定で実験を行うのであれば, ベイジアンでなければならないのである 。

そうではない。いちいち有意水準を調整しない場合も存在する(重回帰分析の各項の有意水準を調整しないのはなぜだろうか?また,相関係数の検定で相関係数の個数を考慮して検定などだれも行わない)。また,どのようなファミリーの検定を想定し,第一種の過誤を統制したいのか,決めていれば多重比較を実行することは問題がないはずである。そもそも,分散分析と多重比較はロジックが異なっており,そこをごっちゃにして議論を行っている。もっといえば,探索的データ解析という方法論も立派に存在している。

第四の利点として,第三の点と 関わるところもあるが,例数設計の問題が挙げられる 。「頑張って 取れるだけデータを取ってみよう 」 という 心理学教育に悪意があるとはいわないが,分析に際して帰無仮説検定を行うのであれば,結論を意図的に作り 出したとの誹りから 免れ
得ないことになる 。

実験であれば,事前に検定力分析をしてサンプルサイズを決定できる。また,多くの場合サンプルサイズが小さく,検定力が不足している研究が多い(脳科学など)。

ベイジアンにとっては得られたデータから 作られる 事後分布は常に一つであり , それが確率的に変わるものではない。 たくさんデータを取ればより確信度が高く , あまりデータが取れなければより 確信度が低くなるだけであるから ,実験者の努力は反映されるし, そのことが結果に影響を与えるものではない。

事後分布はデータの関数で,データが確率的にばらつくのだから,事後分布は確率変数になるはずである。検定においても,データが多くなれば標準誤差が小さくなる(≒日常用語として確信度が高くなる)。「そのことが結果に影響を与えるものではない。 」という結果とは何か?

しかし例数設計を報告し,認可を受けてから 実施するという 作法に則るのが難しい実践場面は少なくないだろう 。 もしなんらかの事情で,事前の例数設計をすることが困難な状況にあるのならば,帰無仮説検定で結論を出すのではなく , ベイジアンでなければならない。

たしかに,仮説検定を行ううえでは,例数設計を行うことが望ましいが,それができなくても,集まったサンプルのサイズから検出できる効果量の大きさは算出できる。その効果量の大きさが事前に想定していたそれと比較することも可能だと考えられる。回帰分析の回帰係数の検定や尤度比検定を使う場合には例数設計を行わないことのほうが多い。また,「ベイジアンでなければならない」という言い回しがよくわからないが,例数設計ができないことがベイズ統計を強要させる根拠にはならない。

第五に, ベイジアンは p値の解釈のように仮想世界の非現実的な値について一喜一憂するのではなく ,事後分布から 自分の仮説が正しい確率を読み取ることができる点である 。解釈が直観的であり ,誤解を生じにくい点は教育的効果が非常に大きいだろう。

帰無仮説検定は基本的に保守的な(慎重な)立場をとっているがそのあたりも考慮はされていなさそうである。事後分布も様々な仮定のもと”非現実的な事後分布”を考えているようにも思う。仮想とか,非現実とか,何を意味しているのだろうか。直感的であることを礼賛しているが,「事後分布から 自分の仮説が正しい確率を読み取ることができる点」ということについては,非常に注意が必要であることを見落としている。あるパラメタの事後分布が特定の区間に入る確率というのは,「真のパラメタがその区間に入る確率」を意味しない。データ生成分布が特定できるコイン投げのような状況であれば,表が出る真の確率を考えてもいいのかもしれないが,一般には,あくまで,「仮想的なモデル,事前分布,データを与えたもとでの事後分布」に過ぎず,それがtruthを反映しているかどうかは事後分布だけからは全くわからない(だから汎化誤差の評価やモデルの評価が必要になる)。そのため,

直観的な発想に沿った解釈ができるということは,誤解を生みにくいことにもつながり ,誤用や悪用を防ぐことにもつながる 。

という主張は非常に危ういものであると考えざるをえない。誤用・悪用のオンパレードになりうる(少なくとも私は人を騙すことができると思う)。


本稿でいうベイジアンとは「ベイズ主義者」 といった人の思想信条に帰属されるものではなく ,考え方の違い, あえて言えば考え方の作法についての話なのである 。

ベイズ統計学は事前分布とデータ生成モデルを考えたもとで算出される予測分布を使って真の分布を推測しよう,という推測方法であり,こういった主義や考え方の違いということを強調するのは望ましくないと感じられる。渡辺先生も述べているように,「正しい推測方法は存在しない」し「最尤推定ベイズ推定も誤っている」のである。

Sober(2008) は統計学における 三つの立場を,異なる問いに立脚したものとして分けて考える 。頻度主義的統計学,特に帰無仮説検定の考え方は,差があるといってよいかどうか,効果があると 認めてよいかどうかという「どちらの道を進むべきか」 についての意思決定の問題である 。 ベイズ主義の考え方はうえで見たように,「どの程度信じられるか」 という 信ぴょう 性の問題である 。第三に尤度主義という 考え方があり , これは手元のデータにもとづいて「採択するべき 仮説はどれか」 という 選択の問題である 。 それぞれ問いの形が違い,必要とする答えが異なるのだから , そもそも 直接比較できるものではないともいえる 。

今までのベイズ推しはなんだったのかと脱力した。このような考え方の違いを述べてどうするのだろうか。結局主義だけでは先に進まないと思うのだが。。。

これはもちろんベイズの公式であるが, モデル, データ , パラメータが全て含まれており , データを取ってモデルやパラメータを考えるという意味で,一般的な科学的営みを表現していると 言えるのではないだろうか。

ベイズの公式をみて,科学の一般的な営みを想像するのは結構だが,最尤推定だって,モデル,データ,パラメタがすべて出てくるし,最尤推定値は所与のモデルのもとで真の分布に(カルバックライブラ距離の意味で)一番近いパラメタを選んでいるという意味できちんと推測をしようとしているし,これが一般的な科学的な営みでないわけがない。

また式1について,事前分布を置かない(あるいは(pθ|M)=1とする )

事前分布が一様,あるいは定数に比例すると書いたほうがよい。

その意味で, ベイズ主義的モデル式は頻度主義の考え方を包括的に含んでいると 言える

包括的だからといって,つねに一般的な方がよいわけではない。一般的なものを使いたいなら,一般化線形混合潜在変数モデルなんかをつかうことも考えられるが,それを使いこなせる人はほぼいない。

さて 尤度主義の考え 方はどうであろうか。尤度主義は,差がないとするモデル(M1) と 差があるとするモデル(M2) のどちらがデータに対して当てはまりが良いかを考えるものである 。 これは式(1) でいうところの分母, すなわち (p y|M1) と (p y|M2) の比較を意味し, この両者の比のことを一般にベイズファクターという 。 より具体的で詳しい説明は,岡田(2018) を参照してもらいたいが, ここではこれらの項が先の式(1) に含まれていること , つまりモデル表現の中に尤度主義に必要な情報が包含されていることを 確認してもらえば十分である 。

尤度比検定・・・orz。そして,周辺尤度が式に含まれているからだからどうだというのだろう。わからない。モデル選択に周辺尤度を使います,といえばそれでいいのではないだろうか。

このように, ベイズの公式は頻度主義,尤度主義の考え方が含まれている , より 一般的な表現なのである 。

ベイズの公式自体は,条件付き確率の定義から当たり前に出てくるもので,とくに意味はないと思う。たとえば,「すべての統計的定理は集合論のことばで記述できる。ゆえに集合論の記述がベイジアンの表現より一般的な表現なのである。」ということも正しく。しかし,そんな言及にはなんの意味もないだろう。

ただし,特にベイズ「主義」 として強調する 場合, その特徴はどこにあるかといえば, モデルを尤度と 事前分布の積で表現していること , またデータ 生成過程をモデリングすることから , まだ見ぬデータの予測についてまで言及できるところにある 。 データを取ってその差があるかないかを 判断するということは,顕現化した 目の前のデータに対する 「説明」 に止まっており , それ以上に言及するのはこれまで「考察」 とされていたものである 。

ごく単純な回帰分析だって,最尤推定値をプラグインして予測分布を作ることができる。というか,これが一般的な最尤推定と呼ばれる方法。なぜ,ベイズのみが予測を扱っていると考えているのだろうか。「まだ見ぬデータの予測」については,その分野の人がただ単にやってないから見たことがないだけでは?と思う。

頻度主義を超えて尤度主義に進んでも , それはあくまでも 手に入れたデータ(y) についての比較であり , それ以上は言及するものではない。 ベイズ主義的であるというのは, このモデルから 言える 次のデータについての「予測」 にまで責任をもつことができる 。

尤度主義が頻度主義の上位互換のような表現だが,そういうものなのだろうか。また,周辺尤度を使って「構造」を当てようとするこころみは,汎化誤差を最小化して予測を最もよくするモデルを選ぶこころみと同じではない。予測が最善のように行っているが,ベイズ統計学においてもベイズ自由エネルギー(=周辺尤度)を使って構造を推測しようとする方法がある(WBICがまさにそれ)。全般的に議論がとっちらかっていて読みにくい印象がある。

真にベイズ統計学的であるということは, ただ推定方法の革新や,分布などの概念の導入にとどまるものではない。 モデルによる 現象の説明はもちろん, 現象の予測にまでその判断基準が設けられ, その上でモデルの母体となる 理論の妥当性を求めることになるだろう。少なくとも 技術上は可能なのだから , 後は我々心理学者が予測をしたかったのかどうかをしっかり 問い直し,決意しなければならない。

真にベイズ統計学的であることは,まず主義を主張しないことからスタートするべきでしょう。「その上でモデルの母体となる 理論の妥当性を求めることになるだろう」といっているが,これはただ単に科学的な営み(とはなんなのか私も明確な答えがあるわけではないが・・・)を述べているに過ぎず,ベイズ統計学がどうこうという水準の話ではない。議論の水準が整っていない。そして,「後は我々心理学者が予測をしたかったのかどうかをしっかり 問い直し,決意しなければならない。」といっているが,何を決意するのだろうか。

こうした批判に対して, Kruschke(2014) は上述の「24回投げて 7回は表が出た」 という 例が釘の話だったらどうか, という 。釘が頭の部分を 底にしてピンと 立ち 上がったら 表, そうでなかったら 裏だと 設定し, 24施行中7回ピンと 立ち 上がったとしても , 帰無仮説検定のア
プローチではこの釘が公平な(立ち 上がるかどうかが50% であるという ) 帰無仮説を棄却できない。常識的にはなんらかの仕掛けを疑うところだが,事前になんら 思い込みを持たないというのであれば, そしてデータ 生成過程に言及しないのであれば,公平な釘と 結論するしかない。 このとき 物理的な法則を加味して, そもそも 表が出にくいだろうという 事前分布を 置けば(例えばBeta(θ | 2, 20) などの事前分布を置けば),正しい判断に寄与することができるにも 関わらず。

これの何が問題なのか理解できないのは私の問題のなのだろうか。たんに「公正な」という言葉を誤用しているにすぎない。この場合,「立ち上がる確率が50%」という帰無仮説が棄却できないということは,なんらかのし掛けがあると考えて差し支えないだろう。なぜ,帰無仮説を棄却できないからといって釘に細工がしていないと考えることを否定されないといけないのだろう。あるいは,検定すべき帰無仮説を別の仮説にしたり,信頼区間を計算すれば「整合的な帰無仮説区間」を求めることができる。そして,このときの「正しい判断」とはなんだろうか。

それでも 前分布を恣意的に置くことに関しては, 批判が残るかもしれない。

それ以前にデータ生成モデルがただしいかどうかを考えるべき。

帰無仮説検定のアプローチの利点は,対象がコインであれ,釘であれ, なんらかの心理学的事象であれ,厳格な手続きに沿って検証を行えば誰にでも 結論が出せることである 。 だからこそ 我々心理学者は, 統計学の専門家になるまえに, いちユーザーとして機械的手続きの習得に努めれば良いのであった。 ただこのことが,手続きマニュアルを精読せず,意図せず間違った使い方をしたり意図的に悪用したりするという 悪習慣を招いてしまっている 。自戒を込めていうが,心理学者のこれまでの統計に対する 態度は, マニュアルを熟読せず, わかったつもりで, あるいはどこかで理解を諦めつつ,使えればよいという 不真面目なものではなかったか。 もっとも , この気軽な関係は心理学という 学問の発展には大きく 寄与したはずであり , その意味でも 従来の統計法を一概に批判できないという 弱みを我々は有しているのかもしれない。

帰無仮説検定が使えないようであれば,ベイズ統計的な方法はもっと使えないのでないか。結局何をみたら,分析が終わるのか,その基準も主体的に決めていかなければならないようなベイズ統計学は自由にできる程度は高いかもしれないが,それがどれくらい分野に寄与するのかはわからない。また,データ収集法(調査方法,実験方法)と解析を区別しているが,学者の専門性というのはその2つをきちんと結びつけるところにもあるのではないかと思う。統計をきちんとやっていないから,クライシスが起こったのだろうか。適当な研究計画やブラックボックスプロトコルに起因する部分はないのだろうか。統計を勉強するのは大事だと思うがそれと同じくらい専門的なデータ収集技術や背景理論を充実させることも重要だと思う。

とはいえ, いちユーザーとしての立場以上に,科学コミュニティに関わる 人間のより 一般的なマナーとして,前提と 自らのアイデアを明示するという 振る 舞いは身につけておくべきであろう 。帰無仮説検定を行うのが良い悪いということ 以上に, どういった前提条件にたって, どのような分析手法をとるのかを明らかにする , マナーを身につけなければならないのである 。

マナーというか,論文にきちんと記述するのが普通ではないでしょうか。

探索的な研究の場合は,「どの程度効果が見積もれるのか」 という 程度問題から 始めるのが妥当であり , いきなり 真偽の決着をつける 頻度主義的意思決定問題にするべきではない。 ベイジアンは厳格な例数設計がなくとも , モデルとデータの関係で見積もりを 立てることができる 。 モデル(尤度 × 事前分布)は自由に立案でき ,数式で表現されることで一般性を保って議論することができる 。 また「今日の事後分布は明日の事前分布」(Lindley, 1972) と 言われるように,知識を積み重ねていくことがベイジアンの利点の一つでもあるから , まずやってみるというのであれば, ベイジアンとの相性は特に良いと 言えるだろう 。 このとき ,結果はどのようなものであれ,程度で表現される 。 すぐに効果が「ある /ない」 という 結論ではなく ,大きくある ,小さくある , という 結果の示され方は,努力が無駄になったという 誤解を 避けることができ , ひいては file drawer問題の解決に貢献することになるだろう 。

「探索的な研究の場合は,「どの程度効果が見積もれるのか」 という 程度問題から 始めるのが妥当であり , いきなり 真偽の決着をつける 頻度主義的意思決定問題にするべきではない」といっているが,”いまこのデータから効果があるといえるか”という行動指針を与える意味で検定には意味があると思います。「 モデル(尤度 × 事前分布)は自由に立案でき ,数式で表現されることで一般性を保って議論することができる 。」といっていますが,ベイジアンでなくとも尤度の設定は自由ですし,厳密に考えるならどの分析手法だって数式で書けます。さらに,推測統計を使わない自由だってあると思います。探索的ならなおさら。「このとき ,結果はどのようなものであれ,程度で表現される 。 」といいますが,それは検定でもそうですし,最尤推定を使ったときの標準誤差をみれば程度はわかります。効果量やその信頼区間を見ても同様です。「大きくある ,小さくある , という 結果の示され方は,努力が無駄になったという 誤解を 避けることができ , ひいては file drawer問題の解決に貢献することになるだろう 。」といっていますが,結局小さい効果というのは見たらわかりますし,そのようなものに興味がある人はどれくらいいるのでしょうか。

議論が深まり ,複数のモデルを比較して雌雄を決する必要になった段階まで進むと ,帰無仮説検定の真価が発揮される 。

そうは思わない。むしろ,データが増え,検討すべき仮説が複雑になっていた場合にモデルが特異であっても汎化誤差が見積もれるベイズ統計学の知識が役に立つと考えられる。また,一般に複雑なモデルの推定は難しい。最適化の計算が失敗することもしばしばある。あるいは,計算時間が現実的でないこともある。そのようななかで,答えが出せる方法は役に立つ。ベイズ統計学の推測方法は最尤推定法よりも安定して推定ができることもある。これは推測の精度の良さとは全く別の話である。とはいえ,そうした方法を採用して暫定的にでも答えをだすことは科学的態度として適切だと考えられる。もちろん,データが大きい場合漸近論が成立していると考えられることもあり,場合によっては最尤推定を使ってもよいと思う。場合によってはベイズ推定法のほうが推定に時間がかかることもある。用は選択の問題であり,「この段階にはこの統計学を使うことがよい」という考え方には賛同できない。

再現性の問題も 含めて, これらのマナーを身につけるためには,心理学教育のカリキュラム全体が関わる 必要がある 。 ベイズ統計学は,初学者に対する 教育的側面についても 利点が多い。 2.2節で述べた 例数設計や解釈の容易さに加え, データ 生成メカニズムを考える 自発性と創造性を育むことができる 。 また「なぜそのような確率分布を用い, その検定統計量に変換しなければならないのか」 といった疑問を封殺し,計算手続きだけを教え込むことによって,嫌悪感を 増幅させるようなこともない。帰無仮説検定の緻密なロジックの組み立ては,針の穴一つ開いた 途端に瓦解するような 危険性を 持っており ,初学者にとってはより 間違いにくい利用法としてのベイズ統計の方が適していると 考えられる 。

しかし,特異なモデルや識別できないモデルについて考えないモデルを考えてしまうこともあるだろう。そういった点を考慮していないのでないか。「また「なぜそのような確率分布を用い, その検定統計量に変換しなければならないのか」 といった疑問を封殺し,計算手続きだけを教え込むことによって,嫌悪感を 増幅させるようなこともない。」といっているが,精密標本論自体が本来高度なので,数学的背景がないと詳細を追うことができない。個人的には,”MCMC”という妖精さんの仕事で,なーんとなく事後分布がもとまるというのも相当に疑問だと思う。また,事後分布を求めるためのプログラムを習うことに抵抗はないのだろうか。「帰無仮説検定の緻密なロジックの組み立ては,針の穴一つ開いた 途端に瓦解するような 危険性を 持っており」といっているが,そもそも帰無仮説検定はケース・バイ・ケースで必要な検定を作ってきたので特殊な場合にしか使えないのは当たり前といえば当たり前である。その前提を犯してしまうとだめなのは,教育すればよいようにも思う(教育してもわかってもらえないのは知っている)。

もっとも , ベイズ統計教育を広めるためには,教える側にも 結果を 読み取ったり 方針を 示したりできるような,技量や能力が求められる 。

これも,結局,先行研究で使われているモデルの焼き直しや,公開されているプログラムをちょっと変える程度なので,これまでの仮説検定とやっていることは変わらないようにも思う。

またベイズ統計学の初学者からよく 聞くこととして,独創的なモデルを作るためにはどのように学べばよいのか, というものがある 。

目的は独創的なモデルを作ることなのだろうか。メカニズムを明らかにしたい現象が先にあるのではないだろうか。なんのためのデータ解析なのか,モデルをいじって遊ぶためではあるまい。

ベイズに限らずモデリングという 観点からテキストを探す必要があるのが現状である 。

テキストはいくらでもあるように思われる。

この ROPEの幅をどれぐらいに置くのかについては,各研究領域で議論して定めねばならない問題である 。 ベイジアンとて, アプリオリにこの大きさであればよいと 決めてかかるわけにはいかないのである

ベイジアンなんだから好きに決めたらいいんじゃないかな。なぜここだけ分野の様子を見てるのだろう。そして,実質的な差の話が今後の方向性なら,効果量でいいんじゃないかな。簡単で解釈しやすいし。


まとめとして気になる点

  1. ベイジアンベイズ統計という言葉をきちんと定義しないまま議論を進めている。
  2. 同様に,頻度論的アプローチなるものが何を意味しているのかも定義がなされていない。
  3. わかりやすい=誤用が少ないと考えている
  4. ベイズでない方法は仮説検定しかないと考え(特に実験研究しか念頭においていない印象),検定を使わずに予測分布について考えている人たちや,行列を分解している人たちの存在を無視しすぎている。
  5. 検定論は確かに欠点はあるが,それでもデータから意思決定をする上で有益な情報を提供していることを軽視している。
  6. 全般的に記述に疑問が多い
  7. 結局,タイトルにあるようなことは述べていない。

MplusのMixture modelでのfactor scoreのout putについて

F列は全体のfactor score,C_Fはmost likelyなクラスに対するfactor score。
あと,fscoreはEAP

http://www.statmodel.com/discussion/messages/13/741.html


Rick Sawatzky posted on Tuesday, June 03, 2008 - 10:10 am
I am fitting an polytomous unidimensional IRT mixture model (using MPlus 5.1 - as in ex.7.27 but with oridinal variables) and requested savedate = cprob. I have three questions:
(1) Is the "F" column in the saved data file the expected a posteriori factor score?
(2) What is the difference between the "F" and "C_F" columns?
(3) Is it possible to save the maximum likelihood scores?
Thank you!

Linda K. Muthen posted on Tuesday, June 03, 2008 - 11:07 am
1. Yes.
2. C-F is for most likely class membership.
3. No.

お金関係のサイト

お金関係で便利なサイト

http://disclosure.edinet-fsa.go.jp/

http://tskl.toushin.or.jp/FdsWeb/view/FDST000000.seam

http://csr-toshokan.net/index.php

http://www.jsda.or.jp/

http://www.jsda.or.jp/shiru/index.html

冊子わかりやすい
http://www.jsda.or.jp/manabu/publications/files/startbook.pdf

https://www.shiruporuto.jp/public/document/container/stat/

https://www.shiruporuto.jp/public/data/movie/hyakka/

https://www.shiruporuto.jp/public/document/container/navi/content/reading264.html
https://www.shiruporuto.jp/public/document/container/navi/content/reading273.html


https://www.toushin.or.jp/investmenttrust/about/what/

https://www.fsa.go.jp/ordinary/index.html#yakudatu

ETFについて
https://www.nikkoam.com/products/etf/about

インデックス投資の本について
https://investment-by-index-invest.com/basic-investment-book-step2/#


年末調整
https://www.nta.go.jp/users/gensen/nencho/index.htm

確定申告
http://www.nta.go.jp/taxes/shiraberu/shinkoku/kakutei.htm
https://www.freee.co.jp/kb/kb-kakuteishinkoku/beginner/

Rの出力の色などを変更したい+forなどの繰り返し出力を1行で表示したい

GitHubにある,crayonパッケージを使えば簡単


github.com


詳しくはヘルプを見たり,exampleを走らせればすぐにわかる。
簡単に使えるので,良さそう。

  library(crayon)
  cat(blue("Hello", "world!\n"))

# Crayon defines the %+% string concatenation operator, to make it easy to assemble stings with different styles.

  cat("... to highlight the " %+% red("search term") %+%
      " in a block of text\n")
# Styles can be combined using the $ operator:

  cat(yellow$bgMagenta$bold('Hello world!\n'))
# See also combine_styles().

# Styles can also be nested, and then inner style takes precedence:

  cat(green(
    'I am a green line ' %+%
    blue$underline$bold('with a blue substring') %+%
    ' that becomes green again!\n'
  ))
It is easy to define your own themes:

  error <- red $ bold
  warn <- magenta $ underline
  note <- cyan
  cat(error("Error: subscript out of bounds!\n"))
  cat(warn("Warning: shorter argument was recycled.\n"))
  cat(note("Note: no such directory.\n"))

# See Also
# make_style() for using the 256 ANSI colors.

# Examples

cat(blue("Hello", "world!"))

cat("... to highlight the " %+% red("search term") %+%
    " in a block of text")

cat(yellow$bgMagenta$bold('Hello world!'))

cat(green(
 'I am a green line ' %+%
 blue$underline$bold('with a blue substring') %+%
 ' that becomes green again!'
))

error <- red $ bold
warn <- magenta $ underline
note <- cyan
cat(error("Error: subscript out of bounds!\n"))
cat(warn("Warning: shorter argument was recycled.\n"))
cat(note("Note: no such directory.\n"))


下のコードを実行すると,iの繰り返しに伴った主力が1行で更新されながら表示される。
(通常は,繰り返しごとに1行出力する。)
ポイントは,flush.console()と,cat関数の頭にある"\r" これでOK。
反復アルゴリズムの途中経過の出力をしたい場合には,よい。

for(i in 1:10) {
  cat("\r i = ",i, ": ",format(Sys.time(),"%H:%M:%S"))
  Sys.sleep(0.5)
  flush.console() 
}

Word(for mac)ショートカットキーメモ

ショートカットキーをいじるといろいろと都合が良い。
→ よく使うやつは,登録しておくこと。慣れてきたら増やすこと。
Windowsのやつとの整合性をつけること。覚えることが多くなるのは嫌だけど。。。)

Word ショートカットメモ
command + cont + 1 黒 (Color: で選択
command + cont + 2 赤
command + cont + 3 青
command + cont + A 網掛け (DefaultCharShadin

すべてのコマンドから
command + cont + E ハイライト(Highlight
command + cont + = 数式 (EquationInsert

フォント
command + option + T Times New roman

デフォルト
option + 左右 一単語分移動