図解・ベイズ統計「超」入門 あいまいなデータから未来を予測する技術 (サイエンス・アイ新書)
ベイズの入門書。
amazonのレビューにあるように非常にわかりやすい。読む価値のある本だ。
この本の良さは、通常とは少し違った説明をして理解を深めてくれることだと思うので、ベイズ統計を学ぶ際に最初に読むのではなく、数式も少しは入った本を読んだうえで、読むのがよいように思った。特に乗法定理を面積で説明するところは非常に分かりやすい。
知っていれば特に問題ではないが、モンティ・ホール問題の回答を1/3(これは模範解)だけ提示するのは理解の躓き石になりそうだ。事前分布を一様分布にした場合、最頻値は1/2であり、直観解である1/2が間違いであるわけではない(1/3は平均値)。本の説明順序や内容から、説明できることには限界があるが、模範解だけの提示はベイズ的ではない気がした。
Mplusで潜在クラスは分析できたが、その後どうするかという問題がある。
例えば、計算した潜在クラスを他のモデルの独立変数や従属変数に使用したいといった場合だ。
最近は、同時推定(共変量を伴った潜在クラス分析やDistal Outcomesを伴った潜在クラス分析やその拡張モデル)を使用するので、ここで示しているのは簡易法である。主成分分析の結果を別の分析に使用するのも同じような戦略である。
話がそれたので潜在クラス分析に戻ろう。
例はUCLAのidreから。
LATENT CLASS ANALYSIS | MPLUS DATA ANALYSIS EXAMPLES
https://stats.idre.ucla.edu/mplus/dae/latent-class-analysis/
データはこちら→ https://stats.idre.ucla.edu/wp-content/uploads/2016/02/lca1.dat
TITLE: Fictitous Latent Class Analysis. DATA: File = https://stats.idre.ucla.edu/wp-content/uploads/2016/02/lca1.dat; VARIAVBLE: NAMES = id item1 item2 item3 item4 item5 item6 item7 item8 item9; USEVARIABLES = item1 item2 item3 item4 item5 item6 item7 item8 item9; CATEGOLICAL = item1 item2 item3 item4 item5 item6 item7 item8 item9; CLASSES = c(3); ANALYSIS: TYPE=mixture; PLOT: type = plot3; series = item1 (1) item2 (2) item3 (3) item4 (4) item5 (5) item6 (6) item7 (7) item8 (8) item9 (9); SAVEDATA: file = lca1_save.txt ; save = cprob; format = free; OUTPUT: tech11 tech14;
DATA
File:ローカルの場所を指定することが多いが、ウェブ上のものも参照できるようだ。
Variable
name:
変数の名前をつける。データとして指定されているdatファイルには変数名がついていない。Mplusでは行頭に変数名を入れれず、変数名はスクリプト内で指定する。usevariables:
分析に使用する変数を指定する。categorical:
カテゴリカル変数を指定(この指定がない場合は潜在プロフィル分析となる)。classes = c(3): クラス数の指定
Analysis
混合モデルを指定。
他の値はデフォルト値で設定されている。
デフォルトは少ないので最終の結果等では指定した方がよい。例:カッコ内説明
STARTS = 200 10(ランダム初期値。デフォルトは10 2)
STITERATION = 50(初期最適化の反復回数。デフォルトは10)
Plot
各項目への条件付き応答確率をプロットで出力する。
type = plot3;を指定する。
seriesオプションを設定し、変数名(プロットの順番)[item1 (1)]と書き、プロットに必要な変数すはすべて指定する。
以下のようなプロットが出力でき、クラスの解釈に有用である。
このプロットは自動的に出力されない。
メニューから[GRAPH]->[View Graphs]->[Estimated probabilities]と選択すると出てくる。
SAVEDATA
各個人の各クラスへの所属確率の推定値と所属クラス番号をファイルを出力できる。主成分分析をした時にデータが出力されるのと同じようなものである。file:
出力するファイル名(場所とファイル名を両方指定する。"C:\Mplus\output.txt"というような感じである。)
save:
cprobとはprobilityesのこと。cはおそらくclassの略なので、クラス確率(Class Probility)ということだろう。cprobは潜在クラス分析でのみ使うコマンドである。.txt(テキスト形式)や.sav(SPSS形式)、.dat形式やcsv形式で出力できる。ただ、出力されるファイル内容は同じなので一工夫必要である。SPSSは多形式取り込みの手順が必要である。csv形式のファイルをSPSSに取り込むという手順と同じだ。
他の形式はエディタで加工する。スペースが半角スペースが5つデータ間に挟まった形で出力されるので、半角スペース5つをカンマに置換するとよい。"_____"→","
このように処理した上で、拡張子をcsvに変更するとExcelなどで開けるようになる。
変数名を加えた出力はファイルはこちら。
ここで出力されたデータを他の分析に使用することができる。
Mplusの説明書に掲載されている例7.7から。
https://www.statmodel.com/usersguide/chapter7.shtml
潜在クラス分析に使用する変数が順序づけられていないカテゴリカル変数である例。
パスダイアグラムは以下のように書く。潜在クラス数が複数あってもCは一つだけである。
TITLE: This is an example of a LCA with unordered categorical latent class indicators using automatic starting values with random starts DATA: FILE IS ex7.7.dat; VARIABLE: NAMES ARE u1-u4; CLASSES = c (2); NOMINAL = u1-u4; ANALYSIS: TYPE = MIXTURE; OUTPUT: TECH11 TECH14;
クラス数の指定のオプション。ここでは2クラスでの計算を指定している。
名目変数として扱う変数を指定する。
カテゴリカル変数の場合にはCATEGORICALと指定する。この分析(順序づけられていないカテゴリカル変数)なのでNOMINALが指定されているが、CATEGORICALと書いても分析結果は基本的には変わらない。実際に走らせてみると、クラス１とクラス２が入れ替わったが、それ以外の出力は同じ。
ただし、CATEGORICALと書かないと反応確率が出力されない。
デフォルトの標準誤差はロバスト推定が使われている。いわゆるロバスト最尤推定法(MLR)である。他の推定法を使用する際には、ESTIMATORオプションを追加する。例えば、ベイズ推定の場合には"ESTIMATOR = BAYES;"を追加する。といっても、MLR以外の推定法として指定できるのは最尤法(ML)とベイズ推定法(BAYES)のみである。ロバスト重み付き最小二乗法(WLSMV)などの最小二乗法系の推定法は使用できない。潜在クラス分析は最尤推定法で行うので当たり前といえば当たり前である。
各オプションは下記の４つ。
ランダム初期値: STARTS = 200 10(デフォルトは10 2)
初期最適化の反復回数: STITERATION = 50(デフォルトは10)
ブートストラップの反復回数: LRBBOOTSTRAP = 50;
ブートストラップの尤度比検定の初期値: LRTSTARTS = 0 0 100 20;
デフォルトに含まれていない出力を指定するコマンド。
Mplusの例ではTECH1とTECH8が指定されているが、ここではTECH11とTECH14を指定している。
TECH11はVuong-Lo-Mendell-Rubin LIikelihood Ratio Test (VLMR) の出力、TECH14はBootstrapped Likelihood Ratio Test (BLRT)の出力ができる。これらの統計量はクラス数が妥当か否かを判断する際には必要となる。どちらか一つ選ぶのであればBLRTなので、TECH14をだけいれることが多いかもしれない。
TECH14によって指定したクラス数より一つ少ないカイ二乗値の差と検定が出力される。P値が低い場合(例えば5%水準以下)であれば、クラス数の増加によってモデルが最適化されていると判断する。下記の出力結果をみると、クラス数を2に指定して、P値は0.000であるので、一つクラスを増やして、CLASSES = c (3);として再分析してみた方がよいかもしれない。TECH14を指定するとブートストラップを施行するので計算時間が割と長くなる。
ちなみに、TECH8オプションはモデルの最適化の履歴が出力されるほか、分析中に所要時間がわかる。ただ、出力が長くなり、使わない部分が増えるので不要かもしれない。
Information Criteria Akaike (AIC) 42850.964 Bayesian (BIC) 42961.756 Sample-Size Adjusted BIC 42907.736 (n* = (n + 2) / 24) FINAL CLASS COUNTS AND PROPORTIONS FOR THE LATENT CLASS PATTERNS BASED ON ESTIMATED POSTERIOR PROBABILITIES Latent Classes 1 3713.04544 0.74261 2 1286.95456 0.25739 CLASSIFICATION QUALITY Entropy 0.343 TECHNICAL 11 OUTPUT VUONG-LO-MENDELL-RUBIN LIKELIHOOD RATIO TEST FOR 1 (H0) VERSUS 2 CLASSES H0 Loglikelihood Value -21445.523 2 Times the Loglikelihood Difference 74.083 Difference in the Number of Parameters 9 Mean 14.753 Standard Deviation 9.570 P-Value 0.0008 LO-MENDELL-RUBIN ADJUSTED LRT TEST Value 73.129 P-Value 0.0008 TECHNICAL 14 OUTPUT PARAMETRIC BOOTSTRAPPED LIKELIHOOD RATIO TEST FOR 1 (H0) VERSUS 2 CLASSES H0 Loglikelihood Value -21445.523 2 Times the Loglikelihood Difference 74.083 Difference in the Number of Parameters 9 Approximate P-Value 0.0000 Successful Bootstrap Draws 10
最も有名な一致係数である。扱うことができるのは名義尺度・順序尺度である。評価者が2名のものが基本だが、κ係数は3名以上に拡張されたバージョンも存在する。Cohenのκ係数というと評価者が2名の場合を指す。
κ係数の目安としてよく使われるのはLandis and Koch (1977)によって示されたものだろうか。
0.0〜0.2: わずかに一致(slight agreement)
0.21〜0.40 まずまずの一致(fair agreement)
0.41〜0.60 中等度の一致(moderate agreement)
0.61〜0.80 かなりの一致(substantial agreement)
0.81〜1.0 ほぼ完全、完全一致(almost perfect or perfect agreement)
Krippendorff (1980)にも基準が掲載されている。
0.67未満 評価しない(discounted)
0.67〜0.80 不確かな結果(conclusions tentatively)
0.80以上 明確な結果(definite conclusions)
Landis and Koch (1977) に比べて厳しい内容である。ともあれ、一致度は0.8以上が目安になるという内容である。
もちろん、この基準には根拠が明確にあるわけではないが、0.8という基準はよく使われている。また、論文などに記載するときには、一般的に認められている水準というような「これが私たちの常識だよ」的な記載がよく見られるが、これらの文献を引用した方が恰好がつくように思う。
CohenのκはSPSSで計算可能である。
評価者を列、ケースを行のデータセットを作る。
[分析]→[クロス集計表]で、図のように行列を設定し、統計量の中の「カッパ」をチェックすると出力される。
なお、SPSSで計算できるのはノーマルなCohenのκで、重みづけκは出せない。
CohenのκはRのirrパッケージで計算ができる。irrパーッケージは重みづけκも計算ができるので必要がある場合は、Rを使うことになるだろう。
irrパッケージ
https://cran.r-project.org/web/packages/irr/index.html
https://cran.r-project.org/web/packages/irr/irr.pdf (pdf)
library(irr) # irrの読み込み
rater01<-c(1,2,1,2,1,1,2,2,1,2,1,1,1,2,1,2)
rater02<-c(2,2,2,1,1,1,2,1,1,1,2,1,2,2,2,2)
d1<-cbind(rater01, rater02)
kappa2(d1)
重みづけのオプションは"equal", "squared"の2つ。"unweighted"がデフォルト値で重みづけなしを意味している。
kappa2(d1, "squared")
equalはすべてのレベルで評価者間の不一致が均等に重み付けされ、squaredは不一致は完全一致からの自乗距離に応じて重み付けされる。
重みづけκの使用を検討するのは、評価が多段階である場合だ。
日本語で以下に解説がある。
http://www.med.osaka-u.ac.jp/pub/kid/clinicaljournalclub12.html
以下、数式を記載する。
κは次のように定義される。
ここでのは観察者の間で相対的に観察された一致度であり（正確度と同一）、
は観察されたデータを用いて各カテゴリーをランダムに観察者の確率を計算する仮説の確率一致である。
は、カテゴリーを
, 項目の数を
と
、評価者の数を
、予想されるカテゴリーを
とすると下記の数式となる。
文献
一致率を計測する際に使ったクリッペンドルフのα係数(Krippendorff's alpha coefficient)のメモ。
一致率というとCohenのκ係数が有名だが、κ係数は評価者が2名という限定がある。
クリッペンドルフのα係数は1)評価者が3名以上に対応し、2)比例・間隔・順序・名義尺度それぞれのオプションが用意されているという点で優れている。
Krippendorff's alpha
https://en.wikipedia.org/wiki/Krippendorff's_alpha
日本語の解説はないが英語は割と充実している。
他の一致率の統計量と同じで、-1〜1の間の値をとり、1が完全一致で、-1が完全不一致である。
Krippendorffのαを出すにはいくつかの方法が用意されている。
・SPSS, SAS, Mplus
http://afhayes.com/spss-sas-and-mplus-macros-and-code.html
Andrew F. Hayesによってマクロが作成され無料で公開されている。
Download KALPHA: kalpha.zip をダウンロードする。zipファイルの中に3つの統計パッケージのマクロが入っている。
・R
パッケージirrで計算ができる。
https://cran.r-project.org/web/packages/irr/index.html
https://cran.r-project.org/web/packages/irr/irr.pdf (pdf)
基本形は以下のようなコマンドラインだ。
kripp.alpha(x, method=c("nominal","ordinal","interval","ratio"))
分析例
# irrの読み込み
library(irr)# データの作成
d1<-c(1,1,1,1,2,1,1)
d2<-c(5,4,5,5,4,4,5)
d3<-c(4,5,4,4,4,4,5)
d4<-c(1,2,1,1,1,1,1)
d5<-c(1,1,1,1,1,1,1)
d6<-c(4,3,3,4,2,4,1)
d7<-c(4,4,4,4,4,4,4)
d8<-c(1,1,2,1,1,1,1)
d9<-c(1,1,1,1,1,1,1)
d10<-c(3,2,3,3,3,3,3)
d11<-c(4,3,4,2,4,4,1)
d12<-c(1,1,1,1,1,1,1)
d<-cbind(d1,d2,d3,d4,d5,d6,d7,d8,d8,d9,d10,d11,d12)# Krippendorffのα係数の算出(順序尺度として計算)
kripp.alpha(d,"ordinal")Krippendorff's alpha
Subjects = 13
Raters = 7
alpha = 0.838
評価者は行にして、ケース列としてデータを入力する。分析例だと7人の評価者でケースは12例ということになる。
α係数は0.838とそこそこ高い値が出ている。
そこそこ高いと書いたものの、どの数字以上が高い数字なのかは、特に基準があるわけではない。研究分野によってだいたいの相場感があるので、それを目安にすることになるだろう。
分析例は順序尺度として計算をしたがオプションを"nominal"にすると名義尺度の計算ができる。比例尺度、間隔尺度の場合も同様にオプションを変えれば計算が可能である。
野村総一郎 メディア用語としての"新型うつ病"のその後 (特集 臨床現場から見た精神疾患の変貌) 臨床精神医学 45(1), 37-42, 2016-01
中島は日本人に備わっていた「徳」（善や正義をつらぬく人格的能力）と「仁」（おもいやり， いつくしみ）が失われたことを指摘する。幼少児期には誰でも自己愛的なものだが，大人に成長してそれを克服し，徳や仁，忍耐を身に付ける。しかし，そのプロセスが現代人に欠如し自己愛が肥大したこと，自分の外からストレスがやってくる，という外在化の考え方が強くなり，自助能力が低くなったこと，などを指摘し，総体として日本人の精神力が弱まったことを原因の1つ（後の2つは精神病理学の衰退と，SSRIを中心とした薬物療法中心主義）としている。これは非常にストレートな言い方であるが，「よく言ってくれた！」と膝を打つ人も多いのではなかろうか？
中島の「日本人の精神力が弱まった」に野村総一郎は「よく言ってくれた！」と膝を打ったようだが、これでは、100時間働いて身体の調子が悪くなった人を「心が弱い」と言っているのと同じである。あまりにもひどい言いようである。
そもそも日本人の精神力が弱まったエビデンスでもあるのだろうか。
昔の日本人には徳や仁があったというのは、近年の右翼的な言動である。徳が「善や正義をつらぬく人格的能力」であり、その能力が日本人にあったならば、昔の日本人は東アジアの国々を相手に戦争をして、人殺しなどしなかったと思うのだが。それとも東アジア人はおもいやり、いつくしむ対象ではなかったのであろか。
パス図を描くソフトを探していると、いろいろな手段があることがわかった。
その一つ、Pencilを紹介。Nightly buildだが、Windows以外にも対応しているようだ。
PENCIL PROJECT http://pencil.evolus.vn/
試しにパス図を書いてみた。
わりと簡単に作成できる。パーツを移動すると、矢印も一緒についてくるなど操作性もいい。
出力形式はpngとsvgの２種類だ。通常はpngで十分だと思う。svgはベクター画像なので、曲線がカクカクになりにくいなど有利な場合もあるが、Wordなどで使うには少し厄介だ。svgを使う場合には、Inkscapeで一度svgファイルを開いて、オブジェクトとしてコピーしてWordに貼り付けるという手順を踏むできるらしい(参照)。
LaTexの場合もInkscapeが必要だそうだ。もちろんpngで保存していればLatexでそのまま使えるので特に悩む必要はない。
svgの場合は、Inkscapeで開いてから、epsかpdfに形式を変換して使用するそうだ(参照)。であれば、最初からpngでやった方が効率的だろう。
