Rでオッズ比と調整済み残差を出す

2x2のクロス集計表の場合はオッズ比と調整済み標準化残差は似たような統計量になる。どちらが好まれるかは分野によって異なるので、適宜、好まれる方を使用するのがよいと思う。社会学ではそれほど見かけることはないが使用されるはずである。オッズ比はなぜか知らないがあまり好まれない傾向にあるように感じる。ちなみに医学ではオッズ比が好まれ、調整済み標準化残差を見かけたことがない。

オッズ比は95%信頼区間を併記するのが普通である。調整済み標準化残差は5%基準、1%基準、0.1%基準をアスタリスク等で表記する。95%信頼区間はどちらかが1を越えていないかを一つ一つ確認しなければならないので、表を見るのがめんどくさい。一方で、調整済み標準化残差のアスタリスク表記は、ぱっと見で分かるので便利だと思う。

有名な統計パッケージでは、SPSSでは調整済み標準化残差は計算できるが、オッズ比の計算ができない。stataはマクロを入れると計算できたように記憶している。とりあえず、計算が一番簡単で効率よくできるのはRであるのは間違いない。

下準備

クロス集計表のマトリックスを作成しておく。

例1 データマトリックス

d1 <- matrix(c(2,15,30,12), nrow=2)

データセットからも作成しておく。

例2 CPS1985のデータを使う

library(AER) #AERパッケージの読み込み
data(CPS1985) #CPS1985の読み込み
d2 <- table(CPS1985$gender, CPS1985$married)

性別と結婚ステータスを選んだ。2x2のクロス表が作りたいのでどちらも二値データの変数である。実際の分析としてはあまり意味がないかもしれない。

調整済み残差

調整済み標準化残差はRのデフォルト機能で出力できる。

d1res<- chisq.test(d1) #カイ二乗検定
d1res$stdres # 調整済み残差

結果は以下。

   [,1]      [,2]
[1,] -4.166099  4.166099
[2,]  4.166099 -4.166099

絶対値4.17である。

d2res <- chisq.test(d2) #カイ二乗検定
d2res$stdres # 調整済み残差

結果は以下。

    no       yes
male    0.259397 -0.259397
female -0.259397  0.259397

こちらは絶対値0.26である。

調整済み標準化残差は標準正規分布区間推定値と同様の値をとる。つまり、95%信頼区間は-1.96~1.96であり、残差ののように期待値からは外れていることを示す数値はこの95%信頼区間の外でああればよい。したがって、1.96より大きいと5%水準で有意な差があることになる。どうように、1%水準は2.58、0.1%水準は3.29、0.01%水準は3.89以上になる。

したがって、例1は4.17であったので0.01%水準で有意だったが、例2は1.96に満たないため、有意ではないということになる。

標準正規分布(余談)

ちなみに1.96というような数字は標準化しているために使用できる基準である。標準正規分布表からこの数字は読み取れる。

正規分布の上側確率
http://aoki2.si.gunma-u.ac.jp/CGI-BIN/tgxp.html

リンク先のでは1.96が例に上がっている。1.96のセルは0.0250であるで、これを2倍(両側にするため)すると0.05=5%水準となる。

オッズ比

オッズ比を出すには、パッケージのインストールが必要である。いくつか存在しているが、一番いろいろできるのはepitoolsパッケージである。

library(epitools) #epitoolsパッケージの読み込み

mid-p法(median-unbiased estimation)

いわゆるオッズ比というとこの方法のことである。

oddsratio.midp(d1) #midpをつけた形で
oddsratio(d2) #oddsratioのみで

単純にオッズ比といったときにはmid-p法のことなので、oddsratio()で構わない。 d1の方を出力すると以下のようになる。

$data
          Outcome
Predictor  Disease1 Disease2 Total
  Exposed1        2       30    32
  Exposed2       15       12    27
  Total          17       42    59

$measure
          odds ratio with 95% C.I.
Predictor    estimate      lower     upper
  Exposed1 1.00000000         NA        NA
  Exposed2 0.05953432 0.00773934 0.2558031

$p.value
          two-sided
Predictor    midp.exact fisher.exact   chi.square
  Exposed1           NA           NA           NA
  Exposed2 3.413834e-05  3.46431e-05 3.098556e-05

$correction
[1] FALSE

attr(,"method")
[1] "median-unbiased estimate & mid-p exact CI"

必要なのは$measureのところなのでそこだけ取り出して出力する。

oddsratio(d1)$measure #後ろに$measureをつける

そうするとオッズ比と95%区間のみが表示される。

          odds ratio with 95% C.I.
Predictor    estimate      lower     upper
  Exposed1 1.00000000         NA        NA
  Exposed2 0.05953432 0.00773934 0.2558031

オッズ比が0.059、95%信頼区間が[0.008, 0.256]である。 クロス表の関連を総当たりで確かめることをよくするが、そういう時に便利だと思う。

その他の方法(参考)

詳しくはこちらを参照。

oddsratio.fisher(d1) # Fisherの方法
oddsratio.wald(d1) # Waldの方法
oddsratio.small(d1) # 正規分布近似+小標本補正

強姦被害者に処女かどうか聞くことは必要か?

警察官が被害者に「処女ですか?」と聞く必要はない - キリンが逆立ちしたピアス」で処女膜損傷が「強姦致傷」と認められた判例のことを知り、判例を調べてみた。

この問題はもともと、山口敬之氏からレイプ被害を訴えた伊藤詩織さんが警察で「処女ですかと聞かれた」という話がこの話の発端であった。この事件から「処女の被害であれば強姦致傷になるからだ」という判例が話題になったようだ。

結論から言うと「処女の被害であれば強姦致傷」は法的には正しい。正確に言えば、処女膜の破損が認められると強姦致傷が認められ刑が重くなる。

もちろん、だからといって処女か否かを聞くのが正しいわけではないし、現在の刑法の判例が正しいということではない。犯罪被害者、警察、法律家の間に認識の差がかなり存在して、警察や法律家にとって処女かどうか聞くことは基本的で重要なことなのだろうが、犯罪被害者にとってはセカンドレイプにもなりうることであり、両者の溝は埋まりそうにない。

問題の「処女膜」は裁判でどの程度扱われているのかを調べてみた。判例データベースLEX/DBで「処女膜」をキーワードにして調べてみたところ128件のヒットがあった。強姦とは無関係の裁判が1件であり、現時点では127件のヒットがある。この数字が多いか少ないかは文脈に依存するが、しばしば裁判において処女膜への言及が確認できるといったところだろうか。

処女膜破損により強姦致傷になった判例

まずは、近年、処女膜破損により強姦致傷になった判例から。処女膜の破損を裁判で扱うのは、古い判例昭和34年10月28日最高裁判所第二小法廷にみられるものという話があったため、状況は今でも変わっていないことを示すためだ。

この事例は処女膜の破損によって強姦致傷になった事例である(津地方裁判所平成27年4月21日, LEX/DB文献番号:25506308)。

同人の両手首をビ ニールテープで椅子の両手すりに縛り付けるなどの暴行を加え,上記一連の暴行・脅迫により,同人の反抗を抑圧し,同人と性交し,その際,同人に全治約1週間を要する処女膜裂傷の傷害を負わせ

罰条に強姦致傷とする刑法181条2項(177条前段)が加わっている。強姦に至る過程で暴力が振るわれていると認定されているので、これらの行為で強姦致傷も成立する可能性が高いが、処女膜裂傷の傷害があるのであれば、強姦致傷は確実に認定されることになる。

成人の事件

処女膜の破裂に関する判例は未成年のものがほとんどだ。しかし、未成年のケースではなく被害者が27歳の事件でも処女膜の破損が証拠として採用され強姦致傷になった事例がある(岡山地方裁判所平成24年9月28日, LEX/DB文献番号: 25483118)。 この判例をみると、処女かどうかの確認は未成年に限定されるわけではなさそうである。伊藤詩織さんも28歳だそうなので、処女膜破損について警察が調べようと考えても不思議はない。

明治時代の判例

さて、処女膜の破損によって強姦が強姦致傷になるという考え方はいつごろから生まれたのかを調べてみた。LEX/DBで最も古い判例である明治44年3月9日大審院の判決(LEX/DB文献番号: 27918395)が以下のものである。

人ノ処女膜ヲ裂傷スルハ即チ人ノ身体ヲ傷害シタルモノニ外ナラス故ニ13歳未満ノ幼女ヲ姦淫スルニ因テ其処女膜ヲ裂傷シタル所為ハ刑法第181条ニ該当スルモノトス

ここで引用されている刑法181条が強制わいせつ等致死傷の条文である。処女膜は身体を傷害したので強姦致傷であるとしている。処女膜破損によって強姦が強姦致傷になるという理論構成はこの時代から大きな変化はないように思える。

処女膜破損は強姦の証拠になるか

少し話題が変わるが、判例の中で「処女膜がやぶれること」が強姦の有力な証拠として扱われていることが気になった。特に、処女膜が破れていないから強姦はなかったという議論が存在していることだ。

いうまでもなく、処女膜は個人差がある器官であり、性交があったとしても損傷があるとは限らない。にもかかわらず、処女膜の有無が証拠として採用されているのである。もちろん、法律家は不確かであっても処女膜に頼らざるを得ない事情はあるのだろう。他に論証できる証拠があれば証拠として不確かな処女膜について取り上げなくてもよいのだが、強姦事件では証拠が十分に揃うケースばかりではない。そのために処女膜に頼らざるを得ないという事情はあるのだろう。

しかし、強姦の認定を左右する証拠として採用されているのは驚きである。実際の判例を見ていこう。

処女膜が破れていなければ強姦はない

被告人は養女が13歳未満であるにもかかわらず暴行・姦淫をしたとされた事件の再審請求審。この事件は有名な事件である。強制わいせつ、強姦で起訴され、有罪判決(懲役12年)の確定判決を言い渡された請求人が再審請求を行った(大阪地方裁判所平成27年2月27日, LEX/DB文献番号: 15505846)。再審が認められている。根拠となったのは以下の事実。

  • 犯罪事実を目撃したとする供述は全て虚偽であるとする目撃者の新供述
  • 強姦当時の病院診療録によると「処女膜は破れていない」と診断されたため、強姦被害はなかった

判決文から直接引用する。■は判例収録の際に伏字になっているところである。

本件再審請求審において検察官が提出した■医院の■の診療録の写し■によれば,同診療録には,■が,平成20年8月29日,■医院を受診し,「処女膜は 破れていない」と診断された旨の記載があることが認められる。再審請求後に検察官がその写しを入手したという同診療録が,なぜ■,同診療録の上記記載内容 自体には,信用性に疑いをさしはさむべき事情は存在しない。これは,請求人から強姦被害を受けていないとする■の新供述を強く裏付けるものといえる。

処女膜が破れていないことだけで強姦を否定したわけではないが、2つの証拠の1つが処女膜の非破損である。

処女膜が破れていないので強姦はない

同事件の再審が認められ無罪が確定した裁判(大阪地方裁判所平成27年10月16日, LEX/DB文献番号: 25541276)。再審請求同様、女膜の状態が証拠として採用されている。この事件では被害を受けたというBが虚偽の供述をしたと後日認めている。

平成20年8月29日時点において,Bの処女膜は破れておらず,この事実自体,被告人によるBへの強姦がなかったことを如実に示すもの

証拠は処女膜だけではなく、被害の供述が翻されていることも考慮すべきだが、「処女膜が破れていないから強姦はない」と明確に述べられている。

刑法の議論では処女膜の非破損は「強姦がなかったことを如実に示すもの」と位置付けられることもあるようだ。

この事件は、最高裁で確定した有罪判決を覆すためにかなりの努力がされたはずで、再審で被告が無罪になり、結果的によい判決であろうが、その証拠に疑念の余地が残る。

memiscパッケージで分析結果を出力

memiscパッケージのmtableコマンドを使うと分析結果などをきれいに出力することができるようだ。

Package ‘memisc’
https://cran.r-project.org/web/packages/memisc/memisc.pdf

事前準備

回帰分析の結果を出すため、サンプルデータから分析を行っておく。以前のエントリーと同様にAERパッケージに入っているCPS1985というデータを使う。

library(AER) #AERパッケージの読み込み
data(CPS1985) #CPS1985の読み込み
lm0 <- lm(wage ~ age, data=CPS1985) #年齢が賃金を決めるモデル
lm1 <- update(lm0, ~. + gender) # 性別をモデルに加える
lm2 <- update(lm1, ~. + education) # 学歴をモデルに加える
lm3 <- update(lm2, ~. + experience) # 仕事の経験年数をモデルに加える

これで4つの回帰分析が実行された。この4つの結果をmtableコマンドでまとめる。

library(memisc) #memiscパッケージの読み込み
mt1234 <- mtable("Model1"=lm0, "Model2"=lm1, "Model3"=lm2,
                 "Model4"=lm3, #モデル名を指定
          summary.stats=c("sigma","R-squared","F","p","BIC","N")
          #出力する統計量を指定
)
print(mt1234) #テキストで出力する。

出力結果は以下のようになる。

==============================================================
                   Model1     Model2     Model3     Model4    
--------------------------------------------------------------
  (Intercept)      6.167***   6.929***  -4.843***  -1.957     
                  (0.723)    (0.720)    (1.244)    (6.835)    
  age              0.078***   0.085***   0.113***  -0.367     
                  (0.019)    (0.018)    (0.017)    (1.120)    
  gender: female             -2.275***  -2.335***  -2.344***  
                             (0.430)    (0.388)    (0.389)    
  education                              0.827***   1.307     
                                        (0.075)    (1.120)    
  experience                                        0.481     
                                                   (1.121)    
--------------------------------------------------------------
  sigma               5.1        4.9        4.5        4.5    
  R-squared           0.0        0.1        0.3        0.3    
  F                  17.2       23.0       59.9       44.9    
  p                   0.0        0.0        0.0        0.0    
  BIC              3264.5     3243.4     3138.2     3144.3    
  N                 534        534        534        534      
==============================================================

オプション

さまざまなオプションが用意されている。

HTMLで出力する

show_html(mt1234)

HTML形式のファイルを書き出すのは下記のコマンド。

write_html(mt1234,file="mt1234.html")

Rのワーキング・ディレクトリに出力される。

TeX形式で出力する。

toLatex(mt1234)

TeXファイルとして書き出すのは以下のコマンド。

write.mtable(mt1234,format="LaTeX", file="mt123.tex")

テキスト形式で出力

テキスト形式でも出力できる。下記のコマンドはタブ区切り。

write.mtable(mt1234,file="mt1234.txt")

CSV形式で出力

CSV形式だとExcelなどの加工が楽になる。colsep=","を入れてコンマ区切りを指定する。

write.mtable(mt1234,file="mt1234.csv",colsep=",")

独立変数をリラベル

独立変数の変数名が英語であったり略語であったりする場合リラベルできる。

mt1234 <- relabel(mt1234,
                 "(Intercept)" = "(切片)", #独立変数のリラベル
                age = "年齢",
               "gender: female" = "性別",
                education = "教育年数",
               experience = "就労経験"
)

独立変数が変更されていることが確認できる。

=========================================================
              Model1     Model2     Model3     Model4    
---------------------------------------------------------
  (切片)      6.167***   6.929***  -4.843***  -1.957     
             (0.723)    (0.720)    (1.244)    (6.835)    
  年齢        0.078***   0.085***   0.113***  -0.367     
             (0.019)    (0.018)    (0.017)    (1.120)    
  性別                  -2.275***  -2.335***  -2.344***  
                        (0.430)    (0.388)    (0.389)    
  教育年数                          0.827***   1.307     
                                   (0.075)    (1.120)    
  就労経験                                     0.481     
                                              (1.121)    
---------------------------------------------------------

統計量の出力を選ぶ

mt1234 <- mtable("Model1"=lm0, "Model2"=lm1, "Model3"=lm2,
                 "Model4"=lm3, #モデル名を指定
          summary.stats=c("sigma","R-squared","F","p","BIC","N")
          #出力する統計量を指定
)

summary.statsの所の指定は適宜変える。指定できるのは下記の統計量。

R-squared
adj. R-squared sigma
F
p
Log-likelihood
Deviance AIC BIC N

有意のアスタリスクを他の記号に変更する。

mt1234 <- mtable("Model1"=lm0, "Model2"=lm1, "Model3"=lm2, "Model4"=lm3,
           signif.symbols=c("a"=.05,
                            "b"=.01,
                            "c"=.001)
)
======================================================
                  Model1   Model2   Model3   Model4   
------------------------------------------------------
  (Intercept)      6.167c   6.929c  -4.843c  -1.957   
                  (0.723)  (0.720)  (1.244)  (6.835)  
  age              0.078c   0.085c   0.113c  -0.367   
                  (0.019)  (0.018)  (0.017)  (1.120)  
  gender: female           -2.275c  -2.335c  -2.344c  
                           (0.430)  (0.388)  (0.389)  
  education                          0.827c   1.307   
                                    (0.075)  (1.120)  
  experience                                  0.481   
                                             (1.121)  
======================================================

解離性同一性障害の盛衰

解離性同一性障害とは一般的な言葉では多重人格に相当する症候群のことである。この診断名は「虚偽記憶(false memory)」の問題を引き起こし、アメリカでは社会問題となった。

虚偽記憶についてはエリザベス・ロフタスのTEDでのスピーチで簡単に知ることができる。
エリザベス・ロフタス: 記憶が語るフィクション | TED Talk | TED.com

虚偽記憶というのは、体験していない記憶のことを言う。アメリカで問題になったのはセラピストが虐待の記憶をセラピーで植え付けるという行為だ。現在の精神的な不調は幼少期の体験、特に親による虐待によって引き起こされるセラピストの一部は考えていた。実際に面接をしてみると、虐待の記憶がない。実際に虐待はなかったからだ。しかし、それは、虐待経験を消し去るために、記憶が抑圧されているからだと解釈する。セラピストは虐待体験を掘り起こそうと様々な手を尽くす。そのうちに、虐待体験が「さもあったかのように」作りだされていくのである。

今回取り上げるのはジョエル・パリスが解離性同一性障害についてまとめた論評。

Joel Paris,
The rise and fall of dissociative identity disorder.
J Nerv Ment Dis. 2012 Dec;200(12):1076-9. doi: 10.1097/NMD.0b013e318275d285.
https://www.ncbi.nlm.nih.gov/pubmed/23197123
PDFへのリンク

解離性同一性障害の作者はフロイトである

解離性同一性障害ジークムント・フロイトの初期の理論で示されたものである(Breuer and Freud, 1893)。彼は、幼児期のトラウマが後々の「ヒステリア(ヒステリー)」を引き起こすと考えていたためであり、実際にセラピーで解離性同一性障害を作り出していたようである。私たちがフロイトの著作で見ることができる解離症状を起こすヒステリーの症例は、フロイトによって作り出した、もしくはフロイトと患者の共作である。

ヒステリーというのは日常用語のヒステリー(女性か泣いたり叫んだり起こったりすること)とは異なる。現在の診断基準で言えば、その症状は解離性障害、身体症状障害の身体化障害に相当する。日常用語のヒステリーに相当するのは「易怒性(irritability)」であり、稀に「かんしゃく(tantrum)」とも表現される。もともとこの言葉は古代ギリシア語の「子宮(ὑστέρα)」から生まれた言葉であり、女性は感情的で理性的に理解できないといった伝統的価値観と女性固有の子宮という器官が結び付けられた問題の多い概念である。にもかかわらず、ヒステリーという言葉は日本語でも未だに使われるし、英語では"hissy"という言葉があり、どの社会でもかなりしぶとく残っている。

フロイトの著作で解離症状を持った症例を見ると、症例の興味深さに惹かれるだろう。それとともに、人間の理解の越えた部分を垣間見ることができ、精神分析学に奥深さを感じるかもしれない。しかし、精神分析が面白いのは、フロイトが患者と共に症例を作り上げていたからである。高い作家性を持ったフロイトという人物が作り上げた症例が興味深いのは当たり前の話であろう。

精神分析のからくりを知ったうえで、娯楽として楽しむ分には問題はない。しかし、多くの精神分析家は、セラピーを娯楽だとは思っていないし、自分の心理臨床でフロイトのやったことを再現しようとするのだ。その結果、解離性同一性障害を作り出すに至り社会問題を引き起こした。

解離性同一性障害が引き起こした社会問題

フロイトの話に感銘をうけた後世の精神分析家が似たような劣化コピーを作り続ける。症例数はとんどんと増加し、ついに、診断基準に掲載されるまでになった。解離性同一性障害精神分析家によるマッチポンプなのだ。

解離性同一性障害の問題は他にもある。パリスは、多くの家族と人々の生活を破壊したと指摘し、多くの精神分析家がフロイトと同様の間違いを起こしたと指摘している。虐待の虚偽記憶を埋め込まれた患者が、親を糾弾し、多く家庭が破壊された。

アメリカでは、解離性同一性障害と虐待の関係について、ずいぶんと認識が改められたようである。カルフォルニアでは父親の虐待体験の虚偽記憶を娘に植え付けたとしてセラピストに対して訴訟が起こり、父親側が勝訴している(Johnston, 1997)。精神分析家が自己満足で患者を増やして、患者の人生や家族を破壊しつづけることに、社会的制裁がようやく下されるようになったのだ。

エビデンス

トラウマと解離を支持するエビデンスはほとんどない(Lynn et al., 2012)。幼児・児童期の虐待によるトラウマと関連が指摘されるのはむしろ他の精神障害であり、その代表例は境界性パーソナリティ障害(BPD)である(Zweig-Frank et al., 1994)。解離性同一性障害がみられるとされても、他の様々な症状(うつ病、不安障害、パーソナリティ障害、薬物使用など)が併存している。

DSMでは有病率が1%と書かれているにも関わらず、実際に報告されている患者数が非常に少なく、DSMの有病率に根拠はない。また、解離性同一性障害を作りだすセラピー以外では見られないという指摘もある(Piper and Merskey, 2004)。

DSMへの掲載

パリスは著作の中で次のように述べている。

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

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

一旦あるカテゴリーが作られると,確固たる科学的な証拠がない限以外すのは極めて困難であった。これによって,ほとんど使われないような診断でマニュアルが溢れかえるという事態を招いた(Rutterand Uher, 2012)。またマニュアルには,ある時期には許容されても,研究によって否定された長らく使用されていない解離性障害のようなカテゴリーも依然として含まれている。 -ジョエル・パリス『現代精神医学を迷路に追い込んだ過剰診断』44ページ

現在では廃止が検討されている(されるべき)診断名の一つとしている。精神分析家によって症例が作り出されたものであるならば、診断基準に掲載するべきではない。パリスは解離性同一性障害DSMから削除することを(たびたび)主張している。

パリスによれば、解離性同一性障害の見方が変化したのは、精神医学のパラダイム変化が関係している。近年は生物学的な捉え方と薬理学的介入が中心になってきているので、精神分析的が退潮したということだ。その社会的な変化により、解離性同一性障害への信ぴょう性も揺らいできたということであろう。

文献

  • Johnston M (1997) Spectral evidence: The Ramona case: Incest, memory, and truth on trial in Napa Valley. Boulder, CO: Westview Press.
  • Lynn SJ, Lilienfeld SO, Merckelbach H, Giesbrecht T, van der Kloet D (2012) Dissociation and dissociative disorders: Challenging conventional wisdom. Curr Dir Psychol Science. 21:48-60.
  • Piper A, Merskey H (2004a) The persistence of folly: A critical examination of dissociative identity disorder. Part I. The excesses of an improbable concept. Can J Psychiatry. 49:592Y600.
  • Zweig-Frank H, Paris J, Guzder J (1994) Psychological risk factors for dissociation in female patients with borderline and non-borderline personality disorders. J Pers Disord. 8:203-209.

うつ病の自記式尺度について潜在クラス分析を行った論文

今回取り上げるのは下記の論文。

Ulbricht CM, Rothschild AJ, Lapane KL. The association between latent depression subtypes and remission after treatment with citalopram: A latent class analysis with distal outcome. J Affect Disord. 2015 Dec 1;188:270-7. doi: 10.1016/j.jad.2015.08.039. Epub 2015 Sep 1. https://www.ncbi.nlm.nih.gov/pubmed/26384013

論文の概要

潜在クラスを使用した論文。 元データはSTAR*Dであり、潜在クラス分析を行っているのはうつ病の自記式尺度であるQIDS-SRである。STAR*Dが実施されたのは10年くらい前のことだが、新しい分析手法であれば、データがやや古くても論文になるようだ。

この論文から読み取れる潜在クラス分析の特徴

精神医学の論文と潜在クラス分析は相性が良い

社会学もそうだが精神医学でもカテゴリカル変数をよく使用する。例えば「うつ/うつではない」といったような「診断」は2値データである。また精神障害の構造化面接(SCID)は、項目の基準を満たしたかどうかで判断する。DSMでは大うつ病では9項目の症状のうち6項目を満たすと診断される。それぞれの症状の項目は2値である。

潜在変数を求める技法として因子分析が伝統的に使われてきたが、カテゴリカル変数は扱えない。カテゴリカル変数をよく使う分野では潜在クラス分析が非常に強力なツールになりうるのだ。

質問項目が決まると再現性が高い

潜在クラス分析の欠点はいくつかあるが、その一つは計測する質問文やその構成によってクラスの分類が大きく変わる可能性があるということである。そういった点から考えると、この論文ではQIDSを使用しており、質問文は固定されている。質問項目が固定していれば再現可能性も高まるはずである。

分析方法

潜在クラス分析の結果

4つのクラスを導き出している。
1. 軽症
2. 中等度
3. 重症+過食
4. 重症+不眠

ただ、BICで判断すると5クラスモデルが支持されたとのことである。ただ男性が9%、女性が11%だったので4クラスモデルを選択したとのことである。BLRT(Bootstrap Likelihood Ratio Test)を行ったとは書かれていない。日本の文献を読むとBLRTを重視する人が多いように感じるが、クラス数の決定はもう少し自由に決定した方がよいのではないだろうか。

BLRTの検定でモデル改善がなかなか示されず、クラスが8つも9つにもなって解釈ができないというのでは、技法の持ち腐れである。理論的な根拠からモデル数を決定するというやり方が許容してもよいのではないだろうか。

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

f:id:iDES:20170626133500p:plain

Table3は他の精神障害と潜在クラスの関係を探ったものである。技法的には多項ロジスティック回帰分析と同じである。参照カテゴリは軽症(Mild)である。

中等度

男性 45歳以上-、PTSD- 、社交恐怖+
女性 45歳以上-、白人-、社交恐怖+

重症+過食

男性 過食症+、社交恐怖+
女性 過食症+、社交恐怖+

上昇+不眠

男性 45歳以上-、全般性不安障害+、PTSD
女性 白人-、全般性不安障害+、PTSD+、社交恐怖+

インプリケーションとしては重症度によって併存症が変化するというところだろうか。特に、重症の方不安障害の併存が高いのが重要そうである。

distal outcomesを伴った潜在クラス分析

タイトルにもふくまれている"distal outcomes"という言葉。これは、潜在変数を独立変数として顕在変数(観測変数)を従属変数とした分析方法の際に使う用語である。従属変数がdistal outcomesに相当する。このブログではまだ分析方法を紹介していないが、潜在変数を独立変数にして解析をするのは割と難しいのである。

従属変数はcitalopramによる寛解である。日本ではレクサプロ(escitalopram)として発売されている薬の光学異性体である(レクサプロはS(L)体のみ)。STAR*Dでは第1段階目に投薬する薬として使用されている。

distal outcomeを伴った潜在クラス分析にもいくつか方法があるがUlbricht et al.(2015)が使用しているのはLanza et al.(2013)である。僕はよく理解できていないが、ベイズ推定を用いてdistal outcomesの分布を導く方法らしい。

Stephanie T. Lanza, Xianming Tan and Bethany C. Bray, Latent Class Analysis With Distal Outcomes: A Flexible Model-Based Approach Struct Equ Modeling. 2013 Jan; 20(1): 1–26. Published online 2013 Jan 29. doi: 10.1080/10705511.2013.742377 https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4240499/

統計パッケージはLanzaらが作っているSASのマクロPROC LCAにDziakらのdisrtal outcomeマクロを利用して計算したのではないかと思われる。

PLOC LCAマクロ
https://methodology.psu.edu/downloads/proclcalta  

LCA Distal Macro
https://methodology.psu.edu/downloads/distal

f:id:iDES:20170626133550p:plain

QIDSによって4つに分けられたクラスの軽症を参照カテゴリにして、残り3つのクラスの投薬による寛解(HRSD≦7)のオッズ比が示されている。 男性では「過食を伴った重症」で効果が弱く、女性では「過食を伴った重症」「不眠を伴った重症」で効果が薄かった。男性の「不眠を伴った重症」も95%信頼区間の端が引っかかっている程度なので、軽症と重症の寛解率は異なるということであろう。軽症の方が抗うつ薬によって寛解に至りやすい。 また、中等度のうつとの差は出なかった。

インプリケーションとしては、次の2点であろうか。

  1. 軽症に比べ、重症は投薬治療で寛解しにくい
  2. 中等度と軽症は投薬治療の効果(寛解率)に差異はない

うつ病の軽症に対しての投薬の効果に関しては論争的であるので、わりと興味深い結果ではある。

Rで回帰分析のモデルを並べる方法

回帰分析のモデル比較をする表をExcelで作らずRでキレイに整形できるという話を聞いたので少し調べてみた。Rのstargazerパッケージ、texregパッケージできれいに整形できるらしい。いずれも、TeX、HTML、テキストの3種類に出力できる能力がある。

下準備として回帰分析の結果を作る。AERパッケージに入っているCPS1985というデータを使う。

library(AER) #パッケージの呼び出し
head(CPS1985) #CPS1985データの行頭のみ表示
?CPS1985 #データについて説明の表示

このデータを使い4つの回帰分析を作り、それぞれlm0~lm3に格納する。

lm0 <- lm(wage ~ age, data=CPS1985) #年齢が賃金を決めるモデル
lm1 <- update(lm0, ~. + gender) # 性別をモデルに加える
lm2 <- update(lm1, ~. + education) # 学歴をモデルに加える
lm3 <- update(lm2, ~. + experience) # 仕事の経験年数をモデルに加える

memiscパッケージ

テキストでまとめる方法。 memiscパッケージの mtable( ) 関数を使う。

library(memisc)
mtable(lm0, lm1, lm2, lm3)
==============================================================
                     lm0        lm1        lm2        lm3     
--------------------------------------------------------------
  (Intercept)      6.167***   6.929***  -4.843***  -1.957     
                  (0.723)    (0.720)    (1.244)    (6.835)    
  age              0.078***   0.085***   0.113***  -0.367     
                  (0.019)    (0.018)    (0.017)    (1.120)    
  gender: female             -2.275***  -2.335***  -2.344***  
                             (0.430)    (0.388)    (0.389)    
  education                              0.827***   1.307     
                                        (0.075)    (1.120)    
  experience                                        0.481     
                                                   (1.121)    
--------------------------------------------------------------
  R-squared            0.0        0.1        0.3        0.3   
  adj. R-squared       0.0        0.1        0.2        0.2   
  sigma                5.1        4.9        4.5        4.5   
  F                   17.2       23.0       59.9       44.9   
  p                    0.0        0.0        0.0        0.0   
  Log-likelihood   -1622.8    -1609.1    -1553.4    -1553.3   
  Deviance         13635.9    12954.1    10514.6    10510.9   
  AIC               3251.6     3226.2     3116.8     3118.6   
  BIC               3264.5     3243.4     3138.2     3144.3   
  N                  534        534        534        534     
==============================================================

きれいにまとめられる。 memiscパッケージは結果をいろいろな形式で出力なできるそうなので、またエントリーを改めたい。

stargazerパッケージ

stargazerパッケージをインストー

install.packages("stargazer")
library(stargazer)

上の回帰分析をまとめると次のようになる。

stargazer(lm0, lm1, lm2, lm3) # TeX形式
stargazer(lm0, lm1, lm2, lm3, type = "html") #HTML形式

デフォルトにもう少し手を入れていこう。

stargazer(lm0, lm1, lm2, lm3, type = "html", 
          title            = "表1 回帰分析のモデル比較", # タイトルを入れる
          covariate.labels = c("年齢", "性別-女性", "教育年数", 
                               "就労経験", "定数"), #独立変数の名前をいれる
          dep.var.labels   = "賃金"
          )

はてなブログは透明の罫線が表示される困った使用なので、画像として張った。下部にHTML出力のリンクをつけたので、生の出力結果がみたい場合はリンク先を見てほしい。

f:id:iDES:20170625155435p:plain HTML出力

texregパッケージ

次に試すのはtexregパッケージ。デフォルトは次のコマンド。

library(texreg)
htmlreg(list(lm0, lm1, lm2, lm3))

タイトルや変数名などを少し加えてみる。

htmlreg(
  list(lm0, lm1, lm2, lm3),
  caption.above = TRUE,
  caption = "表1 回帰分析のモデル比較",
  custom.coef.names = c("定数","年齢", "性別-女性", "教育年数", 
                        "就労経験")
  )

HTML出力

HTMLの出力はtexregパッケージの方がキレイに出力できるようだ。

RのpoLCAパッケージで潜在クラス分析を行う

Rので潜在クラス分析を行う。

Rで潜在クラス分析ができるパッケージは3つある。

他にMclustというパッケージがあるが、つまり潜在プロファイル分析(連続変数)のパッケージで潜在クラス分析には使えない。

ここではpoLCAパッケージを使用する。

データセットgss82を使用する。このデータはpoLCAパッケージに含まれるのでインストールをする必要がある。ちなみに、潜在クラス分析の古典的な解説書であるMcCutcheon(1987: 30)の図3.1で使用されている例である。

gss82の変数は4つである。カッコ内は変数の値である。

  • PURPOSE; 調査の目的について意見(good/depends/waste of time and money )
  • ACCURACY ; 調査の正確性(mostly true/not true)
  • UNDERSTA; 調査の質問を理解していたか (good/fair, poor)
  • COOPERAT; 面接調査に協力的だったか (interested/cooperative/impatient, hostile )

2クラスモデルを走らせてみる。

data(gss82) #データセットをgss82に設定
f <- cbind(PURPOSE,ACCURACY,UNDERSTA,COOPERAT)~1 #fomulaの設定。潜在クラス分析に入れる変数を選ぶ
gss.lc2 <- poLCA(f,gss82,nclass=2,maxiter=3000,nrep=100) #2クラスモデル

maxiterは反復回数である。デフォルトは1000である。収束しない場合にはmaxiterを多少多めにしておいてもよい。nrepはモデル推定をする回数。デフォルトは1である。初期値はprobs.startで与える。

結果は次のようになる

$PURPOSE
            Good Depends Waste of time
class 1:  0.8953  0.0579        0.0468
class 2:  0.2154  0.2066        0.5780

$ACCURACY
          Mostly true Not true
class 1:       0.6367   0.3633
class 2:       0.0297   0.9703

$UNDERSTA
            Good Fair/Poor
class 1:  0.8327    0.1673
class 2:  0.7422    0.2578

$COOPERAT
          Interested Cooperative Impatient
class 1:      0.8840      0.1043    0.0117
class 2:      0.6478      0.2498    0.1024

Estimated class population shares 
 0.8077 0.1923 
 
Predicted class memberships (by modal posterior prob.) 
 0.8136 0.1864 

統計量やフィッティングの指標は以下のように出力される。

AIC(2): 5592.536
BIC(2): 5658.729
G^2(2): 79.33723 (Likelihood ratio/deviance statistic) 
X^2(2): 93.25329 (Chi-square goodness of fit) 

poLCAでBICは出力できるが、BLRT(Bootstrap Likelihood Ratio Test)を出力することはできないようだ。

Latent Class Analysis (Quantitative Applications in the Social Sciences)

Latent Class Analysis (Quantitative Applications in the Social Sciences)