パスカルの生の倦怠とデュルケームのアノミー

パスカルデュルケームが少し似たようなことを言っている気がしたのでエントリを書き始めたのだが、よくよく確認するとパスカル的な考え方をデュルケームは批判的に取り上げ、嘆いているというのが真相だった。

パスカルの生の倦怠

パンセ (中公文庫)

パンセ (中公文庫)

パスカルの『パンセ』から断章139の一部を引用してみる。

人は、いくつかの障害と戦うことによって安息を求める。そして、もしそれらを乗り越えると、安息は、それが生みだす倦怠のために堪えがたくなるので、そこから出て、激動を請い求めなければならなくなる。なぜなら、人は今ある悲惨のことを考えるか、われわれを脅かしている悲惨のことを考えるかのどちらかであるからである。そして、かりにあらゆる方面に対して十分保護されているように見えたところで、倦怠が自分かってに、それが自然に根を張っている心の底から出てきて、その毒で精神を満たさないではおかないであろう。(断章139)

数ヵ月前、一人息子を失い、訴訟や争いごとで打ちひしがれ、つい今朝がたもあんなにくよくよしていたあの男が、今ではもうそんなことを考えていないのは、どうしたわけだろう。驚くことはない。猟犬どもが六時間も前からあんなに猛烈に追いかけている猪が、どこを通るだろうということを見るのですっかりいっぱいになっているのだ。それだけのことでいいのだ。人間というものは、どんなに悲しみで満ちていても、もし人が彼を何か気を紛らすことへの引き込承に成功してくれさえすれば、そのあいだだけは幸福になれるものなのである。また、どんなに幸福だとしても、もし彼が気を紛らされ、倦怠が広がるのを妨げる何かの情念や、楽しみによっていっぱいになっていなければ、やがて悲しくなり、不幸になるだろう。気ばらしなしには、喜びはなく、気を紛らすことがあれば、悲しみはない。(断章139)

断章139は様々なものが含まれているので、この断章のまとめとはならないが、下記の2つの考え方が読み取れる。

  • 生は倦怠である
  • 倦怠の気晴らしによって忘れることができ、幸福になることができる

裕福で自由な生活をしていたパスカルだから捉えられた視点だろうか。生活に追われている当時の大多数の人間にはこういう発想は生まれにくい。

デュルケームアノミー

アノミーというとデュルケームが有名であり『自殺論』で取り上げられたことで有名である。

自殺論 (中公文庫)

自殺論 (中公文庫)

アノミーamonieという言葉は、ギリシア語の秩序や法を意味するnomosにaがついて秩序がない状態を意味するamonosが起源であり、1887年に哲学者のジャン・マリー・ギュイヨー『非宗教の未来』の中で使用したことから、当時のフランスではある程度認知された言葉だったのではないかと考えられている。

ギュイヨーは肯定的な意味でアノミーという言葉を使用しているが、デュルケームは否定的な意味でも使用している。

自殺論でのアノミーは明確に定義づけられたものでもなく、センテンスごとに意味が少し違ってきていたりするのだが、パスカルとの比較で言えば、下記の部分が比較的近いのではないかと思った。

ともかく進歩を、それも可能なかぎり急速な遊歩を協調する説が、一つの信仰箇条となっていった。しかしまた、不安定さの恩恵をたたえるそれらの説とならんで、その説を生む母胎をなした当の状況を一般化するような、新たな説があらわれつつある。それは、生はいとわしいものであるとのべ、生は快楽よりもむしろ苦悩にみちており、人びとをたんに偽りの魅力によって眩惑しているにすぎないとして、生を告発している。混乱が経済界においてもっとも熾烈をきわめているせいもあってか、この説の犠牲者はこの世界にもっとも多い。  実際、工業、商業は、自殺のもっと多い職業にふくまれている。それらは、自由業の水準に匹敵し、しばしはそれをしのいでいる。とりわけ農業にくらべると、前二者の自殺は多い。その理由は、農業が、古くからの規制力がいまなおもっとも強い影響をおよぼしている産業であって、事業熱がいちばん.浸透していないということによる。かつての経済的秩序の一般的な構造はかくのごとくであったにちがいがない、とよく想像させるのはこの農業である。さらにまた、工業における自殺者のうち、雇主を労働者から区別したならば、右の差はもっと大きくなっていたにちがいない。というかは、もっとも手ひどくアノミー状態に冒されているのは、おそらく雇主だからである。金利生活者の法外な自殺率(百万あたり七二〇)は、もっとも富裕な者がもっとも苦悩に苛まれていることを十二分にものがたっている。それは、人びとに従属をしいるものはすべて、アノミー状態のおよぼす効果を経滅してくれるからである。下層階級は、少なくとも上位の諸階級によって制限を課せられている一線があるから、そのことだけからも、かれらの欲望はより限定されている。ところが、その上位にいかなる階級も存在しないような階級のばあい、背後から引きとめてくれる力がはたらかないかぎり、ほとんど必然的に宙に迷いこんでしまうのである。(p.317-9)

デュルケームが批判的に取り上げる発展と共に現れた「生はいとわしいものであるとのべ、生は快楽よりもむしろ苦悩にみちており、人びとをたんに偽りの魅力によって眩惑しているにすぎないとして、生を告発する」説というのは、パスカルの言っていることに近いように思われる。

デュルケームの分析は100年少し前のものだが、現代日本の分析としてさほど違ってはいないだろう。人はある程度の規制がないと幸せに生きていくことができない。最も規制がない金利生活者の自殺者が最も高い。また、雇用主の自殺率の方が高いとも述べている。デュルケームの言葉でいうと「背後から引きとめてくれる力がはたらかないかぎり、ほとんど必然的に宙に迷いこんでしまう」のである。

近代になって普及した「生」の捉え方が、デュルケームよりも何百年も前のパスカルによって記されているのは、パスカルの特殊な環境によって説明ができるものだろうが、それがパスカルが現代的テーマを考える上でも参考になりうるものになっている原因でもあるのだろう。

ソマーズのD[R]

SatataでソマーズのDの出力の方法を以前のエントリで書いたのでRでの出し方も書いておこうと思う。

ryoureadyパッケージ

ryoureadyパッケージを使うのが最も楽なのではないかと思う。 ord.somers.d関数を利用する。データは「データ作成」以降のスクリプトにある。

library(ryouready)
t1 <- table(d3$Science, d2$Math) # 表形式にする
ord.somers.d(t1)

結果は以下。

Somers' d:
    Columns dependent: 0.518
    Rows dependent: 0.345
    Symmetric: 0.414

注意点としては、表形式にしてから投入しなけばいないところである。

  • Columns dependent
    行変数を独立変数、列変数を従属変数としたソマーズのD。
  • Rows dependent
    列変数を独立変数、行変数を従属変数としたソマーズのD。
  • symmetric
    対称ソマーズのD。

InformationValueパッケージ

InformationValueパッケージでもソマーズのDが出力できる。
但し、こちらのパッケージは2値データが0/1という値でないといけない、列変数を独立変数としたソマーズのDしか出力できないという限定的なの機能になっている。

d3b <- d3 # 新しいデータを作成
library(memisc)
d3b$Science <- recode(d3b$Science, 0 <- 1, 1 <- 2) # リコード
d3b$Math <- recode(d3b$Math, 0 <- 1, 1 <- 2) # リコード
library(InformationValue)
somersD(d3b$Science, d2$Math) # 2値データが前
[1] 0.3244561

データ作成

library(mvtnorm)
set.seed(1234) # 乱数シードの固定
N <-1000 # サンプルサイズの指定
Sigma = matrix(c(1,0.54,0.54,1), byrow=TRUE, ncol=2) # c()の中に相関係数を入れる。
mu <- c(10, 10) # 平均値
d1 <- data.frame(rmvnorm(n=N, mean=mu, sigma=Sigma)) # 乱数の作成
d1 <- d1*5
d1 <- round(d1) # まるめ(四捨五入で整数値に)
colnames(d1) <- c("Math","Science") # 列名を挿入
summary(d1)

4分位。順序変数に。

library(OneR)
d2 <- bin(d1, nbins = 4, labels = c(1,2,3,4),  method = "content")
library(dplyr)
d2 <- d2 %>% mutate_if(is.factor, as.ordered)
colnames(d2) <- c("Math.4","Science.4") # 四分位と変数名を変更
d2_num <- d2 %>% mutate_if(is.factor, as.numeric)
colnames(d2_num) <- c("Math.4n","Science.4n") # 四分位の連続と変数名を変更

2分位(中央値)。Factor型に。

d3 <- bin(d1, nbins = 2, labels = c(1,2),  method = "content")
colnames(d3) <- c("Math.2","Science.2") # 二分位と変数名を変更
d3_num <- d3 %>% mutate_if(is.factor, as.numeric)
colnames(d3_num) <- c("Math.2n","Science.2n") # 二分位の連続と変数名を変更

オッズ比のバリエーション

オッズ比にいくつかの計算方法があるということには少し前から気になっていた。といっても僕の勉強する範囲の研究では古典的なオッズ比以外みたことがなく、何のことかよくわからなかった。

Rのepitoolsパッケージでは4種類のオッズ比が出力できる。何が出力されているのか調べてみた。結果から言えば、僕の知らないオッズ比はBiostatistics / Biometricsなどでの発展によって加わったバリエーションのようだ。

研究分野が遠いので、個人的には古典的なオッズ比しか使わないと思う。しかし、Rのepitoolsパッケージの出力するオッズ比は古典的なオッズ比ではないので、統計パッケージを使用する際には注意が必要である。

シミュレーション

先の結論の述べるためにシミュレーションから始めよう。データCPS1985を使用する。

library(AER)
data(CPS1985)
ct1<-table(CPS1985$gender,CPS1985$married)
library(epitools)

epitoolsパッケージのオッズ比を出力する。

oddsratio(ct1)

結果は以下。

$data

          no yes Total
  male   101 188   289
  female  83 162   245
  Total  184 350   534

$measure
        odds ratio with 95% C.I.
         estimate     lower    upper
  male    1.00000        NA       NA
  female  1.04829 0.7324027 1.502553

$p.value
        two-sided
         midp.exact fisher.exact chi.square
  male           NA           NA         NA
  female  0.7967612    0.8550534  0.7953289

いわゆる古典的なオッズ比は手計算では次のように算出する。

(101*162)/(83*188)
[1] 1.048577

epitoolsパッケージのオッズ比のデフォルトでは1.04829だが、手計算では1.048577。値が違うのである。

epitoolsパッケージのoddsratio関数では4つの関数が用意されている。それぞれの値を計算してみた。

MUE: 中央値不偏推定

oddsratio.midp(ct1)$measure

OR = 1.04829

CMLE: 条件付き最尤推定

oddsratio.fisher(ct1)$measure

OR = 1.048508

UMLE: 無条件最尤推定

oddsratio.wald(ct1)$measure

OR = 1.048577

SSAE: 少サンプル調整推定

oddsratio.small(ct1)$measure

OR = 1.030612

古典的なオッズ比は3つ目の無条件最尤推定と同じの値である。oddsratioのデフォルトは中央値不偏推定のようだ。古典的なオッズ比を使用する際にはoddsratio.wald関数を使わないといけない。ここは割と重要なところだと思う。

glm関数でのオッズ比

glm関数でのオッズ比も計算しておこう。

res <- glm(formula=married ~ gender, data=CPS1985, family=binomial)
exp(coef(res))

結果は以下。

(Intercept) genderfemale
   1.861386     1.048577

glmで計算すると無条件最尤推定と同じ値となっている。実用性だけ考えると、オッズ比をRで出すときにはoddsratio.wald関数を使うという点に注意すれば大丈夫そうである。

ちなみに、基礎分析で2変数間のオッズ比と、一般化線形モデルで主力されるオッズ比を比較するのは両方とも最尤推定法を用いて計算していればOKということになる。つまり、基礎分析でオッズ比を出して、その後に多変量解析をしてオッズ比を見る場合、直接比べてよいということだ。

理論的な説明

後回しにしていた理論的説明を少し書いておこう。下記の本に簡単だが説明があった。

Biostatistics for Epidemiology and Public Health Using R

Biostatistics for Epidemiology and Public Health Using R

Chanの本の54頁からのところだが、表2.4はオッズ比に限らず、リスク比などでも推定法が異なることが示されている。

f:id:iDES:20191015164543p:plain

  • CMLE: 条件付き最尤推定 conditional maximum likelihood estimation
  • MUE: 中央値不偏推定 median-unbiased estimation
  • UMLE: 無条件最尤推定 unconditional maximum likelihood estimation
  • SSAE: 少サンプル調整推定 small sample adjustment estimation

One-sampleとかは何かというと、以下の分類らしい。

  • One-sample measurements: 率、リスク、有病率
  • Two-sample measurements: レート比、リスク比、オッズ比

要するに一次のデータか二次のデータかということのようだ。

オッズ比は大きいサンプルでは無条件最尤推定、小さいサンプルでは残り3つの方法が妥当だとされている。

少数サンプルでは分母が0になるときがある。つまりゼロセルがあるということだが、0による除算は不可能であるためJewell's small-sample adjustmentまたは同様の数学的処理を使用するとよい(Rothman, 1976)。

academic.oup.com

f:id:iDES:20191015164625p:plain f:id:iDES:20191015164635p:plain

UMLE 無条件最尤推定

いわゆるオッズ比。 (129) / (72) = 7.714286

CMLE: 条件付き最尤推定

フィッシャーの正確確率検定を元に条件がつけられている。

f:id:iDES:20191015164646p:plain

SSAE: 少サンプル調整推定

観察数または期待数のいずれかが5未満の場合は、少数標本による推定法(Rothman 1976)を用いる。また、Jewell's small-sample adjustmentを使用するようだ。

Nicholas P. Jewell, 1986, On the Bias of Commonly Used Measures of Association for 2 x 2 Tables - JSTOR

f:id:iDES:20191015164714p:plain

それほど複雑な式ではなく、このアプローチであれば確かに分母にゼロが来てもよさそうだ。ただ、オッズ比の値が小さくはならないかとも思う。

Nicholas P. Jewellが2003年に書いた本では別の方法も提案もされている。

Statistics for Epidemiology (Chapman & Hall/CRC Texts in Statistical Science Book 58) (English Edition)

Statistics for Epidemiology (Chapman & Hall/CRC Texts in Statistical Science Book 58) (English Edition)

f:id:iDES:20191015164658p:plain

分子にゼロが来ても対応できそうだし、こちらだと値も多少マシになりそうな気がする。

トリンテリックス(ボルチオキセチン)と高齢者のうつ病

11月に新しい抗うつ薬トリンテリックス(一般名:ボルチオキセチン)が投入されるようだ。効果はSSRI+リフレックス的な効果で、忍容性が高い(継続率70%)ということのようだ(トリンテリックス | kyupinの日記 気が向けば更新)。

新規の抗うつ薬が登場しても、何に使うの?と思うくらい既に多くの抗うつ剤が存在している。今更、新薬がてでるのかと少し冷やかな気持ちになったが、SSRI+リフレックス的という言葉から高齢者のうつ(late-life depression: LLD)に向いているのではと思いついたので、調べると既に有効性を示す論文が書かれていた。

www.ncbi.nlm.nih.gov

Lu AA21004とされているのがトリンテリックスのことである。薬剤名が決まる前から高齢者うつ病のRCTを行っていたようだ。事前マーケティングの時点で高齢者は主たるターゲットに入っていたのかもしれない。

この研究は8週間のRCTである。トリンテリックス群(5mg/day)、サインバルタ群(60mg/day)、プラセボ群の3群比較をする設計になっている。評価はHAM-D24とMADRSとCGI-S。患者の平均年齢は70.6歳である。

成績は下記のようになっている。

トリンテリックス サインバルタ プラセボ
反応(HAM-D24) 53.2%** 35.2%*** 35.2%
反応(MADRS) 59.7%*** 70.7%*** 35.9%
寛解(HAM-D24) 29.2%* 29.2%** 19.3%
寛解(MADRS) 33.8%* 46.9%*** 20.7%

*: P < 0.05, **: P < 0.01, ***: P < 0.001 vs. placebo.

治療成績はさすがにサインバルタには及ばない。しかし、高齢者のうつ病で有効なのはサインバルタと、プラセボより少しマシなパキシルとマイルドな効果があるリフレックスだけなので、ライバルになる商品が実はサインバルタしかない。サインバルタには及ばなかったとしても、十分に使える薬である。

中止率はトリンテリックス群5.8%、サインバルタ群9.9%、プラセボ2.8%で、トリンテリックス最も忍容性が高いことがわかった。

総合すると、トリンテリックスはサインバルタに効果では劣るものの、忍容性が高く薬に弱い高齢者でも飲みやすい。

今後、高齢者のうつ病での第一選択がトリンテリックスになる可能性は非常に高い。

日本でも高齢者は増加し、それに比例して高齢者のうつ病も増えていくため、トリンテリックスはその時代にフィットした抗うつ薬として十分な売り上げが出せるのではないかという気がした。最初は冷やかな気持ちだったが、この薬はたぶん売れる部類の薬になるだろう。

マーケティング的に書くと陰謀論的でいやらしい感じがするが、患者にとっても、効かない薬を飲まされるよりも、臨床的な効果が確かめられている薬を飲めるようになるということは、健康を取り戻す人が増えることであり、とても良いことである。

LLDについて直近のレビューは以下のものだと思うので、こちらの論文が参考になるだろう。

www.ncbi.nlm.nih.gov

相関係数の比較[Stata]

Stataでも相関係数の計算をしてみた。こちらのエントリ尺度水準に適した相関係数とシミュレーションをStataで行ったバージョンである。

データはこちらからダウンロード、もしくは下部にRでの作成方法を掲載しているので、そのままRで走らせると、dta形式のデータが作成される。

ピアソンの積率相関係数

cor Math Science

0.5605

ピアソンの積率相関係数相関係数を無理やり出してみた。本来は他の技法で導くべきものである。

. cor Math_4n Science_4n // 順序x順序
0.4942
. cor Math Science_4n // 連続x順序
0.5032
. cor Math_2n Science_2n //  2値x2値
0.3763
. cor Math_4 Science_2 //  順序x2値
0.4301
. cor Math Science_2n // 連続x2値
0.4246

適した方法は以下の方法である。

  • ピアソンの積率相関係数: 連続変数と連続変数
  • ポリコリック相関係数: 順序変数と順序変数
  • ポリシリアル相関係数: 順序変数と連続変数
  • 順位双列相関係数: 2値データと順序変数
  • テトラコリック相関係数: 2値データと2値データ
  • 点双列相関係数: 2値データと連続変数

ポリコリック相関係数

順序変数と順序変数の場合。

Stata13以下の場合

Kolenikovさんが作った下記のパッケージをインストールが必要。
polychoric, by any other 'namelist'

10未満の数を自動的にカテゴリカル変数、10以上を連続変数として捉えるように設計されている。

findit polychoric

こちらからインストール。

/* 順序x順序: ポリコリック相関係数 */
. polychoric Math_4 Science_4

Variables :  Math_4 Science_4
Type :       polychoric
Rho        = .5888424
S.e.       = .02450705
Goodness of fit tests:
Pearson G2 = 4.2263446, Prob( >chi2(8)) = .83614614
LR X2      = 4.2295511, Prob( >chi2(8)) = .83584131

Stata16ではパッケージのインストールなしで同じコマンドで出力が可能。(追記:2019/10/16)

ポリシリアル相関係数

順序変数と連続変数の場合。

. polychoric Science_4 Math

Variables :  Science_4 Math
Type :       polyserial
Rho        = .57925495
S.e.       = .02113936

Stata16ではパッケージのインストールなしで同じコマンドで出力が可能。(追記:2019/10/16)

順位双列相関係数

2値データと順序変数の場合。
まずはソマーズのDから。
Roger Newsonさん(http://www.rogernewsonresources.org.uk/stata.htm)のパッケージをインストールする必要がある。

ssc inst somersd, replace

書き方はsomersdと書いて先に2値の名義尺度、2つ目に順位尺度を書き入れる。

.somersd Science_2 Math_4

------------------------------------------------------------------------------
             |              Jackknife
   Science_2 |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      Math_4 |   .5183179   .0294415    17.60   0.000     .4606136    .5760223
------------------------------------------------------------------------------

Stata16ではパッケージのインストールなしで同じコマンドで出力が可能。(追記:2019/10/16)

2値データと順序変数はポリコリック相関係数でもあるので、polychoricパッケージでも計算が可能である。

. polychoric Science_2 Math_4

Variables :  Science_2 Math_4
Type :       polychoric
Rho        = .59908999
S.e.       = .0300931
Goodness of fit tests:
Pearson G2 = 2.8323497, Prob( >chi2(2)) = .24264038
LR X2      = 2.83805, Prob( >chi2(2)) = .2419498

方法が異なるので値が異なる。 Stata16ではパッケージのインストールなしで同じコマンドで出力が可能。(追記:2019/10/16)

テトラコリック相関係数

2値データと2値データの場合。テトラコリックはポリコリックの2値の場合なので、特に何も考えず、polychoricで計算すればよい。

. polychoric Math_2 Science_2

Variables :  Math_2 Science_2
Type :       polychoric
Rho        = .60875007
S.e.       = .03596437
Goodness of fit tests:
Pearson G2 = ., Prob( >chi2(.)) = .
LR X2      = ., Prob( >chi2(.)) = .

Stata16ではパッケージのインストールなしで同じコマンドで出力が可能。(追記:2019/10/16)

Stata13ではテトラコリック相関係数はpolychoricパッケージを使わず計算できた。ただし2値の方が0/1でないと動作しないので、作成したサンプルデータは1/2になっていたので、リコードをする必要がある。

. recode Math_2 Science_2 (1=0) (2=1) // 上書きのリコード

結果は以下。

. tetrachoric Math_2 Science_2

   Number of obs =    1000
 Tetrachoric rho =       0.6087
       Std error =       0.0360

Test of Ho: Math_2 and Science_2 are independent
 2-sided exact P =       0.0000

値はpolychoricパッケージで出力すると同じである。

双列相関係数

2値データと連続変数の場合。

. polychoric Science_2 Math

Variables :  Science_2 Math
Type :       polyserial
Rho        = .58812229
S.e.       = .02720657

イシリアルはポリシリアルの2値の場合なので、ポリシリアルとして計算されている。 Stata16ではパッケージのインストールなしで同じコマンドで出力が可能。(追記:2019/10/16)

点双列相関係数

点双列相関係数はポリシリアル相関係数の一つのバリエーションと捉え方とピアソンの積率相関係数のバリエーションと捉え方があるようだ。ポリシリアル相関係数の一つのバリエーションと考えた方が実用的だが、ピアソンのバリエーションの方も紹介しておこう。ちなみに、ここで実用的といっているのは、連続変数でピアソンの積率相関係数を行った値に近くなるという意味である。

調べてみたところStataでは2つ方法があるようだ。 John A. Andersonによるパッケージ。Stata8のころからあるパッケージらしい。
https://www.stata.com/stb/stb17/sg20/pbis.hlp

help pbis

2値データは0/1でないと作動しないので、リコードが必要。

recode Math_2 Science_2 (1=0) (2=1) // 0/1に変更
pbis Science_2 Math

結果は以下。

------------------+------------------+------------------+------------------+
Coef.= 0.4244          t= 14.8072        P>|t| = 0.0001        df=    998

もう一つはEffect sizesを利用する方法。オプションにpbcorrをつけると、点双列相関係数が出力される。Stata13では標準機能である。

esize twosample Math, by(Science_2) pbcorr
                                   Obs per group:
---------------------------------------------------------
        Effect Size |   Estimate     [95% Conf. Interval]
--------------------+------------------------------------
   Point-Biserial r |  -.4246202    -.4712231   -.3742119
---------------------------------------------------------

ピアソンの積率相関係数と同じ値になるはずなのだが、pbisパッケージの値は微妙に違っており、Effect sizesを利用した方法は数字はあっているが符号が逆になる。2つとも微妙な仕上がりである。

結果

変数 相関係数 相関係数 ピアソンで計算した場合
連続x連続 ピアソンの積率相関係数 0.5605 -
順序x順序 ポリコリック相関係数 .5888424 0.4942
連続x順序 ポリシリアル相関係数 .57925495 0.5032
2値x順序 順位双列相関係数 .59908999 0.3763
2値x順序 順位双列相関係数(ソマーズのD) .5075835 0.3763
2値x2値 テトラコリック相関係数 .60875007 0.4301
連続x2値 双列相関係数(ポリシリアル) .58812229 0.4246
連続x2値 点双列相関係数(ピアソン) 0.4246 0.4246

データの作成

サンプルデータはRで作成している。こちらで作成したものより、点数の分散を大きくしてある。理由は点数が10点以上離れないと、ポリシリアル相関係数の計算ができないからである。

library(mvtnorm)
set.seed(1234) # 乱数シードの固定
N <-1000 # サンプルサイズの指定
Sigma = matrix(c(1,0.54,0.54,1), byrow=TRUE, ncol=2) # c()の中に相関係数を入れる。
mu <- c(10, 10) # 平均値
d1 <- data.frame(rmvnorm(n=N, mean=mu, sigma=Sigma)) # 乱数の作成
d1 <- d1*5
d1 <- round(d1) # まるめ(四捨五入で整数値に)
colnames(d1) <- c("Math","Science") # 列名を挿入
summary(d1)

4分位。順序変数に。

library(OneR)
d2 <- bin(d1, nbins = 4, labels = c(1,2,3,4),  method = "content")
library(dplyr)
d2 <- d2 %>% mutate_if(is.factor, as.ordered)
colnames(d2) <- c("Math.4","Science.4") # 四分位と変数名を変更
d2_num <- d2 %>% mutate_if(is.factor, as.numeric)
colnames(d2_num) <- c("Math.4n","Science.4n") # 四分位の連続と変数名を変更

2分位(中央値)。Factor型に。

d3 <- bin(d1, nbins = 2, labels = c(1,2),  method = "content")
colnames(d3) <- c("Math.2","Science.2") # 二分位と変数名を変更
d3_num <- d3 %>% mutate_if(is.factor, as.numeric)
colnames(d3_num) <- c("Math.2n","Science.2n") # 二分位の連続と変数名を変更

データの書き出し

data <- cbind(d1, d2, d2_num, d3, d3_num)
library(foreign)
write.dta(data, file = "TestScore.dta")

2値と連続変数の関連を示す指標

ROC曲線下の面積の続きである。

医学分野では連続変数とカテゴリカル変数の関連を表現する指標としてROC曲線下の面積が利用される理由はよく知らないのだが、他にも使える指標があるのではないかと思い、候補を並べてみた。

勉強をしているとROC曲線下の面積を使用する理由もわかってくるのだろうと期待しつつ、エントリをしておこう。

点双列相関係数

以前に入れたエントリで解説している。片方が2値、もう一方で連続変数であることが条件である。

データはROC曲線下の面積のものを引き続き使用している。試しに計算してみる場合には、前のエントリのデータ作成のところを参照してほしい。

y3 <- as.factor(y2)
library(polycor)
polyserial(x2, y3)

出力は0.7594117
言うまでもないことだが数字そのものをAUCと値を比較する意味はない。

連続変数の妥当性にはピアソンの積率相関係数を使用するので、点双列相関係数を利用してもよさそうにも思えるが、使われていないらしい。こちらも要検討である。

rmsパッケージ

rmsパッケージを使って統計量を一括で計算できる。

library(rms)
lrm(y ~ x)$stats

出力は以下。

Obs          Max Deriv    Model L.R.   d.f.         P            C                                          
1.000000e+02 7.124702e-08 4.493138e+01 1.000000e+00 2.040623e-11 8.558000e-01

Dxy          Gamma        Tau-a         R2             
7.116000e-01 7.317976e-01 3.593939e-01 4.825788e-01

Brier        g            gr           gp
1.549921e-01 2.170988e+00 8.766938e+00 3.606908e-01

ソマーズのD

Dxy = 0.7116

Somers' D。Somer's Dではない。アメリカ英語ではさまぁ~ずに近い発音になると思うが、日本語ではソマーズと書くことが多いはず。2つの順序変数間の相関の指標である。

順序変数とあるように、ソマーズのDの別名はrank biserial correlationであるようだ(https://www.stata.com/statalist/archive/2009-09/msg00159.html)。順位のない名義2値データと順位データの組み合わせで使う。

原著の掲載誌はASRである。

Somers, R. H. (1962). "A new asymmetric measure of association for ordinal variables". American Sociological Review. 27 (6).

Nagelkerkeの疑似R二乗

R2 = 0.4825788

ロジスティック回帰をすると出力される疑似R二乗でよく見る指標。ナゲルケルケと日本語表記されることが多い。オランダ生まれでオランダにこの姓が多いようなのてで、オランダ語の発音でたぶん間違いないだろう。

Nico Nagelkerke -Wikipedia
Nagelkerkeさんの写真(少し下の方)

グッドマン=クラスカルのガンマ

Gamma = 0.7317976

順位相関係数の一種である。1から1までの値をとって、統計的独立にあれば0となる。

グッドマン=クラスカルのガンマ(γ)Goodman and Kruskal's gamma
https://bellcurve.jp/statistics/glossary/932.html

ケンドールの順位相関係数

Tau-a = 0.3593939

こちらはわりと有名どころ。-1から1までの値をとって、統計的独立にあれば0となる。

ケンドールの順位相関係数 - ウィキペディア(Wikipedia)

最後の方はなくてもよかったかもしれないが、いろいろやり方はあるということだろう。

ROC曲線のサンプルサイズの推定

ROC曲線のサンプルサイズの推定について。質問を受けたので推定方法を書いておこう。

pROCパッケージのpower.roc.test関数を利用する。

aucは目標とするACUの値。任意の値である。
sig.levelは有意水準。デフォルトでは5%となっていて書かなくてもいいし、変えたいときは別の値にするとよい。powerは検出力、つまり「1-第二種過誤」である。
KappaはKappa統計量ではなく、実験群と対照群のケースの比率のようだ。2とすると対照群の数が2倍となる。

library(pROC)
power.roc.test(auc = 0.73,
               sig.level = 0.05,
               power = 0.95,
               kappa=2)

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

ncases = 28.15205
ncontrols = 56.3041

ncasesが実験群、ncontrolsが対照群である。
Cohenのκのことを考えると、比較的少なめのサンプルサイズでOKのようだ。

power.roc.test関数は他の値の算出もできるらしい。

有意水準の算出

power.roc.test(ncases = 28,
               ncontrols = 56,
               auc = 0.73,
               power = 0.95,
               sig.level=NULL)

sig.level = 0.05112969

AUCの算出

power.roc.test(ncases = 28,
               ncontrols = 56,
               power = 0.95)

auc = 0.7305903