本当は効果が全くないときに標本の効果量が任意の値以上を超える確率

変な話だけど，研究者にとって，実験をして得る標本効果量の期待値を最大化する方法はないけど，標本効果量が任意の値x以上を取る確率を最大化する最も合理的な方法はある。標本サイズを小さくすること。

母標準化平均差が0であるとき，たとえば，4人対4人という２群について標準化平均差を求める。大雑把な計算だけど，このときd > .8となる確率は，15%くらい。こんな感じで以下のような表を作った。

1群の人数	0.2以上	0.5以上	0.8以上
4人	39%	26%	15%
8人	34%	17%	7%
16人	28%	8%	1%
32人	22%	2%	0%
64人	14%	0%	0%
128人	6%	0%	0%
256人	1%	0%	0%

もちろん，人数が少ない方が誤差が大きくなって，任意の値以上の標本効果量が出やすくなる。

…そんなに，「有意差はなかったけど効果量は小あった」といって「実質科学的な」議論をするのがお好みだったら，どんどん標本を小さくしたらお好みに沿う。標本が小さい方がそういう可能性はどんどん上がる。

…なんていっててもしょうがないので，ひとこと，
殊更小標本の効果量の点推定値を当て込んで大げさな議論をするのはやめましょう。
効果量の信頼区間も合わせて報告しましょう。

（「どうしたら大きい効果量が得られるんですか」というご質問を受けたことを思い出して）

クラスタリングあれこれ

local shirnkingに基づくクラスタリングっていう手法があるんだって。

http://math.yorku.ca/~stevenw/pub/sw14.pdf

これをやるRパッケージもある。

https://cran.r-project.org/web/packages/clues/clues.pdf

まあ，よさそう。

このパッケージ面白くて，面白いデータセットがある。（まあよくあるタイプなんだけど）

Broken Ringっていうの。こういう2変量のデータ。

クラスタリングの精度を見たりするにはよさげなデータなんだけど，

人間の目に対してのわかりやすさに反して，よく見ると数字的には難しそうなデータよね。

各クラスターの距離とか，歪みとか，そういう分布の形が…

人間の目（人間のパターン認識能力）にはもちろん，問答無用にこう映る。

よほどのことがないと，人間間で分類に差はでないだろう。

5を指定したk-means法だと，

人間からしたらなんでそうなる？？？みたいな。

ちなみに，5を指定した混合分布モデルだと，

なんでそこだけ！！？？的な。

で，まあそのlocal shrinkingによるクラスタリングだと，完璧とかいう話なんだけど。

ううむ。

たとえばirisなら，

正直人間の目ではよくわからない。1クラス離れているのはわかる。

答えは，

こういう答えをみてから，やっと「俺にはそう見えてたけど？」みたいな顔をするのが普通だ。

人間のパターン認識は別に距離の近さだとか，混ざった分布の尤度だとか，必ずしもそういうものでないかもしれない，みたいな。そんなことを感じられる。

このデータの悪意っていうか，その意匠がね。

ううむ。

R連携アプリ開発関係メモ（1）：Rscriptでbat実行するときに，コマンドプロンプトを表示しない簡単な方法はVBScript

美しいリレーの話。

まずはRのスクリプトファイル（.r）を用意する（a.rとする）。

それでこのスクリプトファイルをRscriptで実行するバッチファイル（b.bat）を用意する。

このバッチファイルを実行するVBScript（c.vbs）を用意する。

vbsは，コマンドプロンプトの表示なしでbatをバックグラウンドで実行させることができる。

HSPといったプログラム言語で以上のスクリプトを自動生成する実行ファイルを用意する（d.exe）。HSPでは，vbsを起動させる命令execなどがある。このd.exeからc.vbsを実行させるとよい。

最初のRスクリプトで.csv，.html，または.pdfで分析結果を出力するようにする（knitr，hwriterなど）。さらにRの命令の中でこれらのファイルを即座に自動に開くことができる。この流れ。

作業の流れとしては，d.exeから，a.r，b.bat，c.vbsをそれぞれ自動生成。それでc.vbsを起動し，コマンドプロンプト表示無しでb.batがRscriptに渡してa.rを開始。a.rが.htmlなどで出力を行い，開く。たとえばブラウザでその.htmlが立ち上がった段階で，a.r，b.bat，c.vbsを削除する。すると，まったくRで計算させたように見せないで，実行ファイルだけからすべて計算してブラウザに出力したように見える。

d.exeにUIをつければ，実質はRなんだけど，ぜんぜんRを話しているようには見えない統計アプリが開発できる（ていうかもうできちゃった）。

というか開発技術的にすごく簡単なので誰でもこれでRでやる作業を実行ファイルでまとめて自動化できるし，R言語とかbatとか知らないひと向けのアプリにできる。

頻度主義的・ベイジアン標本サイズ決定

概論的にはこの論文がいいかも。

Adcock, C. J. (1997). Sample size determination: a review. Journal of the Royal Statistical Society: Series D (The Statistician), 46(2), 261-283.

実験する前に，適切な標本サイズを決定しましょう，というそういう話。

基本的には，検定力分析を行う，ていうのも手なんだけど，平均値とかの精度，というか，たとえば信頼区間の幅や信用区間の幅を満たす標本サイズをもとめましょう，という考え方もある。

今日は，平均値の信頼区間および信用区間の幅を満たす標本サイズを求めるほうについて。

たとえば，母数が平均50，標準偏差10である分布を考える。この分布からnケース抽出して標本を取るとき，その標本から見たときの母平均の信頼区間の幅が10になるnはいくつ
だろう？みたいな。たとえば，[45, 55]とか。

もちろん，これは簡単に求まる。

えっと，いちいち方程式で書いていくのめんどいので，Rの関数書いて計算しちゃえ。

n.length<-function(l,s,level=.95){
n<-ceiling(4*s^2*(qnorm((1+level)/2)^2)/(l^2))
n
}

lは，信頼区間の幅，sは標準偏差。

ていうか，SampleSizeMeansというパッケージがある。

https://cran.r-project.org/web/packages/SampleSizeMeans/SampleSizeMeans.pdf

install.packages("SampleSizeMeans")
library(SampleSizeMeans)

mu.freq(10,1/100)

この関数，ちょっと厄介なのはλ（精度）を入れる。λは1/分散。だから前に自分が書いた関数のほうが（自分には）楽。

16人だと。ちょっとブートストラップ的に見てみる。

m<-as.numeric(0)

for(i in 1:1000){
m[i]<-mean(rnorm(16,50,10))
}

quantile(m,c(.025,.975))

で，これには結構問題ある。

信頼区間の幅，は任意に決められる。
でも，分散が既知，ということはまあありそうじゃない。
だから，結構これまで予想をつけてやってたんだけど，予想は結構外れる。
わからんのだから実験しているのに，決まった値で予想をつけるのはむずい。

これを，もっと柔軟にできないか，ってことで。
ベイズ的に，未知の分散について，まったく知らんってことはないだろう，事前分布としてある程度，分布なら予想がつくだろう，って発想になった。もちろんベイズ信用区間の幅で，ということで。

これは，この論文などで知られる。

Joseph, L., & Belisle, P. (1997). Bayesian sample size determination for normal means and differences between normal means. Journal of the Royal Statistical Society: Series D (The Statistician), 46(2), 209-226.

これで紹介している方法は，SampleSizeMeansというパッケージで全部計算できる。
この論文では，いろいろな基準が比べられているのだけど，以下のような基準がある。

• 平均長基準（ALC）
• 平均重複率基準（ACC）
• 最低結果基準（WOC）

詳しくは上の２つの論文を。

ベイズ的に求めるので，まず，分散（精度）についての事前分布を考える必要がある。
この方法では，精度λ（1/分散）が母数（α，β）のガンマ分布だと考える。

ベイジアンがガンマ分布についていうとき，なぜか頻繁に母数をα，βとかというんだけど，一般にαはk，つまり形状母数。βは比率母数のことね，尺度母数θとは違う。比率母数は1/θ。なのでβ = 1/θとして考えれば混乱しない。ガンマ分布の期待値は，αとβでいうと，α/βになる。

というか，よくある事前分布話だけど，1/分散が従うだろう分布がガンマ分布っていうのはまあいい。ただ，αとかβの値なんて思いもつかんがなw

えっと期待値がα／βで，λ = 1/σ^2になのね。とりあえず，λ = 1/10^2を中心として，こんな分布（α = 2，β=200）を成すと勝手に考える。形状母数2で勝手に固定した。本当は事前分布の標本サイズを決めるんだけど，無情報 = 0にする。

で，SampleSizeMeansパッケージでやると，

mu.alc(10,2,200,0,.95)
mu.acc(10,2,200,0,.95)

大体，25人から31人くらいね。
頻度主義的なやり方より，分布で考えている分，広め。

これ当たり前だけど，めっちゃ狭い精度の分布を事前分布に指定してやると頻度主義的ば求め方と一致していく。

たとえば，α = 2000，β = 200000とか。

mu.alc(10,2000,200000,0,.95)
mu.acc(10,2000,200000,0,.95)

そのメタ分析なんだかな

そういう論文ってどれだけの数あるかわからないのだけど，こんなメタ分析があるって考える。
（潜在変数としての）英語の力と（同じく潜在変数としての）国語（母語？）の力の関係をメタ分析でもとめる，という話にしよう。

基本的に，英語の力（の推定値）は，英語の力を表すとする課題の成績の値それ自身ではない。国語も同じである。英語の力（の推定値）は基本的に英語の力を表すとする課題の成績と相関係数があるはずだ，と考える。英語の力っていうのがあり，それが高ければ英語テストに高い成績を取るだろう，と。

そして，世の中に，英語の力を表すとする課題は2つしかなく（TOEICとTOEFLとか），国語の力を表すとする課題も2つしかないとする。

仮に，以下の様な係数をもつ因子分析モデルがあり，これが真のモデルだと考える。

これは，大体こんな相関係数行列にフィットする構造だ。（頑張って探した）

	英語テスト1	英語テスト2	国語テスト1	国語テスト2
英語テスト1	1.00	0.70	0.42	0.40
英語テスト2	0.70	1.00	0.44	0.42
国語テスト1	0.42	0.44	1.00	0.76
国語テスト2	0.40	0.42	0.76	1.00

さて，＜真の＞英語力と国語力の相関係数は.57くらい。

で，この世界のとき，任意に英語テストのどっちかと国語テストのどっちかをランダムに選んで，その組み合わせの相関係数を報告している研究が25件あるとする。

当たり前すぎるのだけど，これらの研究をメタ分析して，つまり相関係数の重み付け平均を求めたとしても，.57にはならない。
せいぜいがそれぞれの観測変数の相関係数であるr = .40ちょっとに近づく。

実際に数値の例を出す。
その25個の研究の標本サイズは，平均50，標準偏差10くらいだとする。
すべての研究の（両方の）課題の信頼性係数は固定的に.80だとする→汚れた観測値が得られる。
英語テストと国語テストの組み合わせは完全にランダムだとする。
こういう条件下の研究で固定効果のメタ分析（OP）をする。

結果はほい。こんな感じ。

区間で言えば，母相関係数の95%信頼区間は[.26, .37]くらい。
希薄化の修正を行えば，[.35, .45]くらい。な。当たり前だけど。

で，この研究をもって，英語力と国語力の相関は.40くらいだ！なんていったらちょっと違う。解釈がよくない。
あくまでも英語テストと国語テストの相関係数として解釈するのならいい。

任意の観測変数の組み合わせに対してメタ分析をしておいて，潜在変数間の関係として結果を解釈するのは，ちょっとなんだかな。

実際の研究はもっと複雑だ。
同じ潜在変数の指標として使われているものが2つところじゃなくてたくさん混ざっているものもある。それらのモデルはしかも得てして不明だ。処遇の効果，というのも基本的には同じだ。医療系のメタ分析のように薬，または処置，そして成果変数が揃っているわけではない。外国語教育研究では，むしろ成果変数（指導の効果を観るテスト）が同じなんてことはほとんどない。結果は無駄だといわない，ただ，解釈が本当にそれでいいのかは悩みどころ。

これで，こういうのを見たら，メタ分析だから科学的だとか，高いエビデンスを示すとか，そういう話もちょっと一歩引いて考えたくなる。メタ分析が優れた手法だというのは間違いない。

crs<-c(42,42,40,44)
n<-as.numeric(25)
c<-as.numeric(25)
c2<-as.numeric(25)
r<-as.numeric(25)
r2<-as.numeric(25)
l<-paste("Study",1:25)
for(i in 1:25){
n[i]<-round(rnorm(1,50,10),0)
c[i]<-sample(crs,1)
c2[i]<-c[1]*sqrt(.80^2)
r[i]<-cor(mvrnorm(n[i],mu=c(50,50),Sigma=matrix(c(100,c[i],c[i],100),2,2)))[1,2]
r2[i]<-cor(mvrnorm(n[i],mu=c(50,50),Sigma=matrix(c(100,c2[i],c2[i],100),2,2)))[1,2]
}

metacor.OP(r,n,l)
metacor.OP(r2,n,l)

ジョンソンのSU分布

ジョンソンの分布っていうのは，こんな感じ。PDFとかCDFとか。

Johnson's SU-distribution - Wikipedia, the free encyclopedia

この分布は4パラミターある。

1. γ
2. δ
3. ξ
4. λ

このように贅沢な関数なので，やはり表現力は高い。
Rでは，SuppDistsパッケージなどでこれを扱える。

https://cran.r-project.org/web/packages/SuppDists/SuppDists.pdf

たとえば，微妙に歪んだこんな分布を考える（ex-Gaussian分布から乱数生成）。

#psychパッケージとretimesパッケージをつかう
library(psych)
library(retimes)

#乱数を生成
d<-rexgauss(1000,3000,200,1000)

#モーメントを確認
describe(d)

#ヒストグラムで可視化
hist(d,col="lightgray")

#ジョンソン分布へフィットさせてパラミタを得る
param<-JohnsonFit(d)

#パラミタを確認
param

#このパラミタにもとづいて乱数を生成
s<-rJohnson(1000,param)

#この乱数をヒストグラムで可視化
hist(s,col="lightgray")

ここからex-Gaussianのパラミタを得てみる。

timefit(s)

私の場合だと，μ = 2998，σ = 219，τ = 976で，正解がそれぞれ3000, 200, 1000だから，まあいい感じだった。

確率密度関数を求めるには，dJohnsonっていう関数を使う。

x<-seq(0,10000,.01)
plot(x,dJohnson(x,param),type="l",ylab="Probability",ylim=c(0,.001))
lines(x,dexgauss(x,3000,200,1000),type="l",ylab="Probability",col=2)

赤はex-Gaussian分布の関数。

いい感じ。

正答率がチャンスレート以上か調べるベイズ因子

ある被験者になんかの判断課題，32試行やってもらったとしよう。

正答数は18個。なので正答率は18/32 = .56くらい。

この人の判断は，チンパンジー（チャンスレート）より優れているとどれくらいいえるのだろうか。ベイズ因子で考える。

まずは，普通に母比率の検定をRでやってみる。

binom.test(18,32)

p = .60くらいなので棄却できない。チャンスレートではないとはいえない。

この関数は区間推定を出してくれる，[.38, .74]くらいだ。

さて次，ベイズ因子で考える。Moreyのパッケージの場合，

bf<-proportionBF(18,32,.5)
bf
1/bf

比率が.50な方に2，そうじゃない方（≠.50）に0.5くらい。

お，まあ，ギャンブル・レベルだけど比率が.50ってのもありそうじゃないか。
でも，これ当たり前だけど，事前分布の設定によってぜんぜん違うから注意。
デフォルトがいつもいいわけじゃない。

比率が0.5より大きいモデルと小さいモデルについてもこんな感じで。

bf<-proportionBF(18,32,.5,nullInterval=c(0,.5))
bf
bf[2]/bf[1]

bf = 2.8くらい大きいモデルに支持的。チャンスレートよりは高いかもとはいえそうだ。

事後分布からのサンプリングもできる。M-Hだって。

bf<-proportionBF(18,32,.5)
bf
chains<-posterior(bf,iterations=10000)
summary(chains)
plot(chains[,"p"])

こんな感じ。

この分布の2.5%から97.5%の範囲は，大体.40から.70くらいだった。

うむ。