井出草平の研究ノート

共変量を伴った潜在クラス分析

潜在クラス分析の潜在変数を従属変数にして分析を行う場合の技法。 パス・ダイアグラムで描くと次のようになる。

f:id:iDES:20171109122256p:plain

例えば、パーソナリティについて質問をして、潜在クラス分析をすると、パーソナリティが4つに分類できることが分かったとしよう。 これらの各パーソナリティに性別が関連しているかを調べたいとする。

この場合、パーソナリティの質問項目はu_1,u_2,....u_nになる。 潜在クラス分析を行ってできた潜在変数がcである。 そして、性別がxであり、共変量となる。

共変量は連続変数でもカテゴリカル変数でも構わない。

ここではMuthénらのセミナーで使われた例を使う。 http://www.statmodel.com/download/Topic%205.pdf

サンプルはこちらから http://www.statmodel.com/examples/webnotes/webnote4.exe

自己解凍機能付きのzipファイルである。使用するのはzip内のasb.datというデータファイルである。

TITLE:     LCA of 9 ASB items with three covariates
DATA:      FILE = asb.dat;
           FORMAT = 34x 51f2;
VARIABLE:  NAMES = property fight shoplift lt50 gt50 force
           threat injure pot drug soldpot solddrug con auto
           bldg goods gambling dsm1-dsm22 male black hisp
           single divorce dropout college onset f1 f2 f3 age94;
           USEVARIABLES = property fight shoplift lt50 threat
           pot drug con goods age94 male black;
           CLASSES = c(4);
           CATEGORICAL = property-goods;
ANALYSIS:  TYPE = MIXTURE;
MODEL:
           %OVERALL%
           c#3 ON age94 male black;
           %c#1%
           [property$1-goods$1*0];
           %c#2%
           [property$1-goods$1*1];
           %c#3%
           [property$1-goods$1*2];
           %c#4%
           [property$1-goods$1*3];
OUTPUT:    TECH1 TECH8;

MODELの%c#1%以下の行はこのスクリプトを動かすだけであれば不要である。%OVERALL%の2行でOKだ。 上のCATEGORICALオプションでカテゴリカル変数がproperty-goodsと指定されているので、潜在クラス分析に使用する変数はこの範囲だと認識しているためである。

とはいえ、もう少し複雑なモデル構築になるとどの変数を潜在クラス分析に含むかはきちんと指定しないといないので、長いモデルの記述の書式も重要である。 %c#1%はクラス1についてであり、property$1-goods$1は潜在クラス分析に含む変数の指定である。最後にクラス1から0を割り振る。潜在クラス分析で導かれた潜在変数が従属変数になるため、共変量の出力は3つのクラスだけ表示される。文末の数字はどの変数に対しての効果かを明確にするために必要とされる。従って、分析をする際に0から始まる数字を適宜振り分ける必要がある。

どのクラス数が適切かは事前に分析しておくとよい。Bootstrap Likelihood Ratio Test(TECH14)で数理的に最適なクラスを導く(参照)のが一般的だが、理論的に想定されるクラスに分けるであったり、所属ケースが非常に少ないので分析に支障があるためクラスを減らすという場合もある。

分析内容の説明だが、asbというのはAntisocial Behaviorの略で、犯罪行為がデータに入っている。propertyは財産犯であり、潜在クラスに入れる変数は犯罪種別ごとの集計結果である。共変量として入るage94=0は16歳、1は17歳、maleは男性(女性=0) black(黒人=1)である。

共変量の部分の出力以下のようになる。

Categorical Latent Variables

                                                    Two-Tailed
                    Estimate       S.E.  Est./S.E.    P-Value
 C#1        ON
    AGE94             -0.285      0.028    -10.046      0.000
    MALE               2.578      0.151     17.086      0.000
    BLACK              0.158      0.139      1.141      0.254

 C#2        ON
    AGE94              0.069      0.022      3.182      0.001
    MALE               0.187      0.110      1.702      0.089
    BLACK             -0.606      0.139     -4.357      0.000

 C#3        ON
    AGE94             -0.317      0.028    -11.311      0.000
    MALE               1.459      0.101     14.431      0.000
    BLACK              0.999      0.117      8.513      0.000

 Intercepts
    C#1               -1.822      0.174    -10.485      0.000
    C#2               -0.748      0.103     -7.258      0.000
    C#3               -0.324      0.125     -2.600      0.009

右列に有意確率が表示されるので仮説検定はそれを利用する。 潜在クラス分析の部分の出力は以下のようになる。

Chi-Square Test of Model Fit for the Binary and Ordered Categorical
(Ordinal) Outcomes

          Pearson Chi-Square
          Value                           1120.441
          Degrees of Freedom                   472
          P-Value                           0.0000
          Entropy                            0.690

FINAL CLASS COUNTS AND PROPORTIONS FOR THE LATENT CLASS PATTERNS
BASED ON ESTIMATED POSTERIOR PROBABILITIES

    Latent
   Classes

       1        928.39438          0.12673
       2       1499.08710          0.20463
       3       2249.49632          0.30706
       4       2649.02219          0.36159

CLASSIFICATION OF INDIVIDUALS BASED ON THEIR MOST LIKELY LATENT CLASS MEMBERSHIP
Class Counts and Proportions

    Latent
   Classes

       1              920          0.12558
       2             1433          0.19560
       3             2154          0.29402
       4             2819          0.38479

パリス『現代精神医学を迷路に追い込んだ過剰診断』

現代精神医学を迷路に追い込んだ過剰診断 -人生のあらゆる不幸に診断名をつけるDSMの罪

現代精神医学を迷路に追い込んだ過剰診断 -人生のあらゆる不幸に診断名をつけるDSMの罪

DSMと過剰診断について書かれた本。 DSMが過剰診断を生み出すという主張はよく見るが、パリスは不適切な使用によって過剰診断は増加すると主張している。

もし精神医学に関する過剰診断が心配ならば,実践する臨床家に注目すべきで,マニュアルに注目すべきではない。端的に言えば, DSM-5自体は問題ではなく,それに過大な価値を置く我々に問題がある。

これは非常に重要な指摘である。 典型例ではない症例が過剰診断になる場合には何らかの原因があるとパリスは考えている。例えば、うつ病双極性障害ADHDの場合であれば、薬物によって治療が可能である誘惑となり過剰診断となる。

PTSDの場合は、祖国のために命を懸けて戦った軍人への負い目や心理療法家が子ども時代に原因を求める傾向などがあげられる。

原書のタイトルはOverdiagnosis in Psychiatry, How Modern Psychiatry Lost Its Way While Creating a Diagnosis for Almost All of Life’s Misfortunesである。邦訳の「DSMの罪」というタイトルはミスリーディングである。

パリスはDSMは聖書では単なる不完全なマニュアルに過ぎないが代わりなるような有望な診断システムは存在せず、問題はあれど使わざるを得ないという立場である。別の言い方では、DSMはコミュニケーシンの手段であるとも言われている。問題は使用者の方にあるということを何度も言っており、よくあるDSM批判とは一線を画している。

PTSDの章は非常に勉強になった。

死に瀕するような闘いの中で人がPTSDになる一方,戦争帰還兵にはしばしば薬物乱用やうつが見られた。このような薬物乱用やうつは必ずしも兵役と関係しているわけではなく,軍に所属する以前から患者が抱えていた問題だった(Young,1997)。またアメリカ退役軍人援護局で兵役関連のPTSDのための治療を受けていた患者の中には,戦闘の経験がない者もいた(McNally, 2007)。しかしこれらの病院にいる患者に治療を受けさせるのには, PTSDは都合の良い診断名であった。

イスラエルの研究者たちによると,ほとんどの兵士は想像を絶するほど過酷な戦闘を経験した後でもPTSDを発症しないという(Zoharet al., 2009)。

一般の人々の心的外傷では,オーストラリアの消防士の前向き研究(McFarlane,1989)と大規模な一般人口(Breslauet al.,1991)で,PTSDの患者の多くが以前から高レベルの神経症的傾向や心的外傷などの脆弱性を抱えていることが明らかになった。したがって,このようにすでに順応性が不安定な人にとって,不運な出来事は「転換点」となる。

つまり、トラウマを受ける以前から心的な問題を抱えていた人が、外傷体験をきっかけに状態がひどくなったり、場合によって、トラウマ体験以前の状態を不問に付して「外傷体験が現在の自分の状態が悪い根源である」という説明を生み出している可能性がある。パリスの本では「心的外傷はあくまでトリガーであり,PTSDを発症する十分な理由にはならない(McNally, 2003)」と書かれている。

そのうえでパリスは「PTSDと診断したいという原動力の背後には,サイコセラピスト(心理療法家)の存在がある。彼らの多くは,臨床現場で心的外傷を中心的課題として見ている」と指摘している。

心理療法家は、何かしらの体験が精神障害の原因であることを想定しがちである。一般的には、幼児期の体験を原因に説明することを好むが、それは思い込みであることが多い。「「抑圧された記憶」を明らかにしようという流行があった。また保育園の職員を無実の罪で刑務所送りにしたりということもあった(McHugh,2005)」。

もちろんPTSDがないという話ではない。トラウマティックな出来事を体験したからといって、必ずしもその出来事が精神状態の不調の原因とは限らず、現在はあまりにもトラウマへの還元がされすぎているという指摘である。

『図解・ベイズ統計「超」入門』


ベイズの入門書。
amazonのレビューにあるように非常にわかりやすい。読む価値のある本だ。
この本の良さは、通常とは少し違った説明をして理解を深めてくれることだと思うので、ベイズ統計を学ぶ際に最初に読むのではなく、数式も少しは入った本を読んだうえで、読むのがよいように思った。特に乗法定理を面積で説明するところは非常に分かりやすい。


知っていれば特に問題ではないが、モンティ・ホール問題の回答を1/3(これは模範解)だけ提示するのは理解の躓き石になりそうだ。事前分布を一様分布にした場合、最頻値は1/2であり、直観解である1/2が間違いであるわけではない(1/3は平均値)。本の説明順序や内容から、説明できることには限界があるが、模範解だけの提示はベイズ的ではない気がした。

Mplusにおける潜在クラス分析で各ケースのクラス所属確率を出力する

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などで開けるようになる。


変数名を加えた出力はファイルはこちら
ここで出力されたデータを他の分析に使用することができる。

Output
参照 http://d.hatena.ne.jp/iDES/20170601/1496265443

Mplusで行う潜在クラス分析(順序づけられていないカテゴリカル変数)

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;
CLASSES

クラス数の指定のオプション。ここでは2クラスでの計算を指定している。

NOMINAL

名目変数として扱う変数を指定する。
カテゴリカル変数の場合にはCATEGORICALと指定する。この分析(順序づけられていないカテゴリカル変数)なのでNOMINALが指定されているが、CATEGORICALと書いても分析結果は基本的には変わらない。実際に走らせてみると、クラス1とクラス2が入れ替わったが、それ以外の出力は同じ。
ただし、CATEGORICALと書かないと反応確率が出力されない。

ANALYSIS

デフォルトの標準誤差はロバスト推定が使われている。いわゆるロバスト最尤推定法(MLR)である。他の推定法を使用する際には、ESTIMATORオプションを追加する。例えば、ベイズ推定の場合には"ESTIMATOR = BAYES;"を追加する。といっても、MLR以外の推定法として指定できるのは最尤法(ML)とベイズ推定法(BAYES)のみである。ロバスト重み付き最小二乗法(WLSMV)などの最小二乗法系の推定法は使用できない。潜在クラス分析は最尤推定法で行うので当たり前といえば当たり前である。
各オプションは下記の4つ。
ランダム初期値: STARTS = 200 10(デフォルトは10 2)
初期最適化の反復回数: STITERATION = 50(デフォルトは10)
ブートストラップの反復回数: LRBBOOTSTRAP = 50;
ブートストラップの尤度比検定の初期値: LRTSTARTS = 0 0 100 20;

OUTPUT

デフォルトに含まれていない出力を指定するコマンド。
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


以下、数式を記載する。

κは次のように定義される。
\kappa \equiv \frac{p_o - p_e}{1 - p_e} = 1- \frac{1 - p_o}{1 - p_e}, \!

ここでのp_oは観察者の間で相対的に観察された一致度であり(正確度と同一)、p_eは観察されたデータを用いて各カテゴリーをランダムに観察者の確率を計算する仮説の確率一致である。

P_eは、カテゴリーをk, 項目の数をNn_{ki}、評価者の数をi、予想されるカテゴリーをkとすると下記の数式となる。

 p_e = \frac{1}{N^2} \sum_k n_{k1}n_{k2}

文献

クリッペンドルフのα係数

一致率を計測する際に使ったクリッペンドルフのα係数(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"にすると名義尺度の計算ができる。比例尺度、間隔尺度の場合も同様にオプションを変えれば計算が可能である。