井出草平の研究ノート

春日武彦は腹の底で何を考えているか

「あけぼのばし自立研修センター」の件で精神科医である春日武彦氏が注目されている。
この訴訟では、成仁病院の医師が逮捕監禁罪で訴えられている。

f:id:iDES:20191202040857p:plain https://twitter.com/pentaxxx/status/1200740076600250368

「あけぼのばし自立研修センター」はひきこもりの「引き出し屋」である。
本人の許諾なく、連れ去り、施設に閉じ込めるのは、拉致と監禁であり、刑事罰の対象となる。

拉致監禁をしたのは「あけぼのばし自立研修センター」であり、「成仁病院」はセンター経由で強制入院を実施したようだ。
入院の際には、法定要件を満たしておらず、全裸にオムツをはかされて3日間身体拘束され、入院期間は50日に及んだと告訴した男性は述べている。

春日武彦の肩書

現在は「顧問」と名乗っており、制度的には「管理者」のようだ。
https://www.nisseikyo.or.jp/hospital_search/hospital_Info.php?id=1224&dt=20161209143629

ただ、今年の秋くらいまでの肩書は院長であった。
https://www.igaku-shoin.co.jp/paperDetail.do?id=PA03335_03
http://rakukai.com/2019/05/29/getsurei-20190615/

春日武彦氏は院長ではないという声がある*1ようだが、問題の訴訟の時には院長だったのだと思われる。 告訴された医師が春日武彦氏なのか不明である。

春日武彦氏ではなかったとしても、病院所属の医師個人が、院長の許可なしに他の団体との提携をしていることは考えにくい。
本人の関与がなかったとしても、管理的な立場にいる者の責任は問われるべきであろう。

春日武彦精神科医は腹の底で何を考えているか』

同書206頁から「引きともりは「心の病」か」というタイトルで、ひきこもりの話が書かれている。

小手先の解決法を模索してみても効果はない。現実的な対応は、まずは家族の硬直した価値観を変えることから始まるだろう。ただしいきなり変えることは困難だし、どうしてもある程度の時間経過が必要だろう、おそらく年単位の。そしてじっくりと時聞をかけて本人と親、双方が「ああ、もっと別な考え方、別な生き方だっであるんだ」と思えるようになって互いに牽制し合うことから脱却した状態を、「和解」と称することになるだろう。

そういった意味では、引きこもりは病というよりは家族病理を和解へと至らしめるためのプロセスと見倣すべきかもしれない(98|和解という形を、ゴールインとして設定することが実際的であると考える医師)。(pp.201)

ひきこもりの支援論として、変なことを言っているようには思えないし、拉致監禁に関わっている人の感じもこの文章からは感じ取れない。

池上正樹×春日武彦 対談

news.livedoor.com

ジャーナリストの池上正樹氏との対談。
「ゴールを設定せずに子供の話を聞いてあげてほしい」と春日武彦氏は述べている。 先の本でもひきこもりは親子の和解がメインテーマであったが、この対談でも同じスタンスで語られている。
ちなみに、2019年8月27日の記事で当時の肩書は成仁病院院長である。

ひきこもりの拉致監禁に手を貸していたとなると「和解」とは程遠い臨床をしていたことになる。
親の依頼で、本人の許諾なしに、拉致監禁をした後に、親子の和解などできるものだろうか。 拉致された側は一生、親を許さないだろう。

行政主催のイベントの講師に

春日武彦氏は東京を中心に、行政主催のイベントで、ひきこもり関連の講師として頻繁に登壇しているようである。
例えば、今年は文京区や大田区で講演をしていることが確認できる。

言動と行動の不一致が明らかであると、春日武彦が「腹の底で何を考えているか」が気になるところである。
次回の講演では「春日武彦は腹の底で何を考えているか」というテーマで、是非、ひきこもり支援の本音を語っていただきたいものである。

*1:斎藤環さんのツイートのレスを参照

Stata LCA プラグインを使用した潜在クラス分析 導入偏[Stata]

Stata15から標準機能で潜在クラス分析ができるようになっている。

ただ、基本的な機能しかないため、Stataの標準機能の潜在クラス分析は現在のところ実用に耐えられない。そこで、面倒な作業が必要になるが、Lanzaらが作成したStata LCA プラグインを使わざるを得ない。

Stataの標準機能

Stataの標準機能の潜在クラス分析は下記のエントリで解説を行っている。

ides.hatenablog.com

Stataの標準機能にかけているものは2つある。

  1. BLRTが計算できない。 潜在クラス分析を使った論文で求められる統計量であるBLRTがStataの標準機能では計算できない。無料のRでも少し工夫をすれば計算できる(参照)ので、有料であるStataを選ぶのは合理的な選択とは言いづらい。

  2. モデルの拡張ができない。 潜在クラス分析はクラス分け(グルーピング)をしただけでは、興味深い結果にならないことが多い。実際には、他の変数との関連性を検討することが必要になる。しかし、Stataの標準機能ではこのモデルの拡張ができない。

ここの2つの理由から現在のStataの標準機能での潜在クラス分析は勧められない。

そこで、Lanzaらの作ったプラグインを使う必要が出てくるのである。ただ、パッケージではなく、プラグインであるため、導入と運用にはハードルが少し高い。このハードルを越えられれば、MplusやLatent Goldなどのソフトウェアと同程度の分析が可能になる。

Stata LCA プラグインの注意点と他のソフトウェア

Windowsにしか対応していない。

  • Mplus Macにも対応。現在の最新技術まで利用できる。

  • Latent Gold Windowsのみ対応。現在の最新技術まで利用できる。GUIで操作性が非常によい。アカデミックでも30万弱、一般で80万とかなり高額。

  • SAS LCA PROC SAS用にLanzaらが配布しているプラグイン。無料である。ただしSASの方が高額である。SASMacにも対応しているが、LanzaらのプラグインWindowsのみ対応。もちろん無料版のSAS University Editionには対応していない。

Macの人は1) Mplusを使う(他のソフトウェアの使用をあきらめる)、2) デュアルブートWindowsを動かす、3)Winodws機を購入するという3点が選択肢になる。

Stata LCA プラグインのダウンロード

www.methodology.psu.edu

上記のページから32-bit版/64-bit版のファイルをダウンロードする。
Windowsが64bitに対応していれば、どちらでもよいが、BLRTをおなうことを想定すると、計算時にメモリを多く使用するため、メモリの使用に制限がない64-bit版の方がおそらくよいだろう。

プラグインのインストール

ダウンロードしたファイルを解凍

ファイルを常時設置するディレクトリ(フォルダ)が必要になる。ファイルをどこに置くかはよく考えた方がよいだろう。マニュアルの例は次の場所。“D:\project\Stata_lca\Release\” 。
一般的にこのようなソフトウェアは日本語など2バイト文字を含むアドレス上に置くと動作しないことが多いため、Windowsのユーザー名を日本語で入れている人はユーザーフォルダは避けた方がよい。

下記で変更を行う"lca.do"、"doLCA.ado"は解凍フォルダに入っているため、解凍したファイルすべてを自分で決めた場所のフォルダに移動する。

doファイルの設定

doファイルは毎回設定をする必要があるので、後で説明。

"doLCA.ado"ファイル内のパスを変更

  1. ファイル"doLCA.ado"を開く。 注:Stataの[ファイル]→[開く]から。ファイルをダブルクリックしても開かない。
  2. コードの一番最後の行(869行目)の"lca.dll"のパスを先ほどと同じように変更する。 例では“D:\project\Stata_lca\Release\lca.dll”。必ず"lca.dll"までを含めたパスを書く必要がある。
  3. 変更を保存。

doファイルの設定

同梱されている"lca.do"ファイルは5つの例題が掲載されている。動作確認するには例1だけ走らせるのがよいだろう。

  1. "lca.do"ファイルをコピーし、"lca_ex01.do"など任意の名前にリネームする。
  2. 22行目(//matrix list r(rhoSTD))までが例1なので、任意の名前に変更したファイルの23行目以降削除する。これで例1だけのdoファイルになる。
  3. コードの5行目で、パス“D:\project\Stata_lca\Release\”をさきほどファイルを置いた場所に変更する。この行には“/CHANGE THIS PATH TO MATCH THE FILE LOCATION ON YOUR MACHINE/”という文字が書きこまれてあるのですぐに発見できるはずだ。
  4. 変更を保存し、doファイルを実行する。

結果

例1はこれでとりあえず走るはずである。

scalars:
                 r(df) =  4031
        r(EntropyRsqd) =  .8585102062978721
         r(EntropyRaw) =  227.7190382066842
        r(AdjustedBIC) =  979.1363404006836
               r(caic) =  1246.403913372541
                r(bic) =  1182.403913372541
                r(aic) =  868.3075755176842
           r(Gsquared) =  740.3075755176842
      r(loglikelihood) =  -3366.110782824834
          r(iteration) =  591

matrices:
             r(rhoSTD) :  24 x 5
                r(rho) :  24 x 5
           r(gammaSTD) :  1 x 5
              r(gamma) :  1 x 5

. matrix list r(gamma)

r(gamma)[1,5]
       Class1     Class2     Class3     Class4     Class5
r1  .67019096  .06416831  .14781815  .06554844  .05227414

. matrix list r(gammaSTD)

r(gammaSTD)[1,5]
       Class1     Class2     Class3     Class4     Class5
r1  .02334551  .01335266  .02312087  .00856993  .01521007

所属確率が出力されないのは21行目の"matrix list r(rho)"が"//"でミュートされているためである。
実際に分析する場合には、"//"を取り除く。標準化した結果が欲しい場合には、22行目の//matrix list r(rhoSTD)のミュートを外す。

実行に必要になるもの

2回目以降の分析で必要になるものが3ある。

  1. doファイル スクリプトを書き込むのはdoファイルである。

  2. データファイル 同梱の"LcaSampleDataset.txt"を開くとデータの構造がわかるはずである。変数名は入れず、データは数字で入力。データとデータの間は半角スペースで区切る。

以下のような形式である。

1 3 15 2 2 2 2 2 2 2 2 2 2 2 2 1
1 1 17 2 2 2 2 2 2 2 2 2 2 2 2 2
1 0 17 2 2 2 1 2 2 2 2 2 2 2 2 3
2 2 19 2 2 2 1 2 2 2 2 2 2 2 2 4
1 4 17 2 2 2 2 2 2 2 2 2 2 2 2 5
2 2 15 2 2 2 2 2 2 2 2 2 2 2 2 6
2 5 15 2 1 1 1 1 2 2 2 2 2 2 2 7
2 1 11 2 2 2 2 2 2 2 2 2 2 2 2 8
1 2 14 1 2 2 2 1 2 2 2 2 2 2 1 9
1 0 17 2 2 2 2 2 2 2 2 2 2 2 2 10

3. 変数名ファイル データファイルと同じ名前で作成する。拡張子はdct。変数名はこのファイルに格納する。

同梱の"LcaSampleDataset.dct"は下記の内容が書かれている。

dictionary using LcaSampleDataset.txt
{
   gender TVHours MomEduc SmokedBefore13 DailySmoke DroveDrunk
   DrankBefore13 BingeDrink MarijuanaBefore13 CocaineEver
   GlueEver MethEver EcstasyEver SexBefore13 ManyPartners ID
}

表示できないので改行をしているが、実際には改行はしてはいけない。 編集はテキストエディタを利用する。こちらもdoファイルと同様に、分析ごと(正確には新しいデータを使うごと)に設定が必要になる。

テキストエディタ

テキストエディタを使ったことがない人には下記の2つのソフトウェアの導入をお勧めする。

多機能のテキストエディタAtomがよい。

atom.io

ただし、Atomはプログラミングからメモまで多くの機能を持った多機能エディタで設定もそれなりに勉強が必要なので、dctファイルの編集のためだけに導入するのはハードルが高い。

シンプルな機能のみに対応したテキストエディタTeraPadがよいと思う。

tera-net.com

AtomTeraPadともフリーウェアである。
秀丸EmEditor、SublimeText3などでもよい。使い慣れた利用するとよいと思う。

トリンテリックス(ボルチオキセチン)と性機能障害

トリンテリックス(ボルチオキセチン)について一つ気づいていなかったことがある。下記の記事を読んで気付いた。

www.mdmag.com

性機能障害が他の抗うつ薬に比べて低いという補足記述が2018年10月22日にFDAに受理されたという記事である。

研究

元になった研究はJacobsen et al.(2015)である。

www.ncbi.nlm.nih.gov

セレクサ、パキシルジェイゾロフトで治療され反応した参加者の薬を無作為にトリンテリックス(10/ 20mg; n = 225)とレクサプロ(10 20mg; n = 222)に切り替えた。期間は8週間である。
性機能の尺度はCSFQ-14(the Changes in Sexual Functioning Questionnaire Short Form)、抗うつ剤の有効性はMADRS、全般的な評価としてCGIを使用した。
性機能障害はレクサプロに比較して少なかった。また、両方のグループで抗うつ薬の有効性は継続していた。

世界で市販されている性機能障害の出にくい抗うつ剤

  • リフレックス・レメロン(ミルタザピン)
  • トリンテリックス(ボルチオキセチン)
  • ブプロピオン(未承認)
  • アゴメラチン(未承認)
  • ビラゾドン(未承認)

性機能障害が少ない抗うつ薬は海外ではいくつか選択肢がある。しかし日本では2つしかない。

リフレックス・レメロンは眠気によって日中のパフォーマンスが低下することなど、副作用が原因で服薬が難しい人が多い。 日本では、高齢者がターゲットなのではとブログで書いたが、性機能障害もターゲットであるのは間違いなさそうである。

研究者コメント

バージニア大学医学部精神医学および神経行動科学科の主任研究調査員であるアニタ・クレイトン医学博士のコメントが掲載されている。

性機能障害は、SSRIを処方されたときにうつ病と闘う最も一般的で厄介な副作用患者の1人です。これらの厄介な副作用を具体的に調べるために研究を設計しました。うつ病の治療効果を失わずに、性的副作用が潜在的に少ない薬剤に変更することは、うつ病患者にとって重要な選択肢です。

トリンテリックス(ボルチオキセチン)のレビュー

トリンテリックスが11月27日発売ということなので、トリンテリックスのレビュー論文を簡単に紹介しておこう。適応の大うつ病性障害以外に全般性不安症にも使用されている。

トリンテリックスのレビューはなぜかたくさんあって、適切なものを選ぶのが難しい。最近のもので内容が広く浅いものということでSowa-Kućma et al.(2017)を抜粋して紹介しよう。薬理の部分ほとんど訳出していない。

www.ncbi.nlm.nih.gov

以下はレビューの抜粋である。

うつ病性障害(MDD)

トリンテリックスは多くの系統的レビューがある(Citrome et al. 2014,Baldwin et al. 2016, Berhan et al. 2014, Pae et al. 2015, Meeker et al. 2015, Fu et al. 2015, Thase et al. 2016)。統合された結果によるとトリンテリックスによる抑うつ症状の改善はあまり大きくない(約0.2の標準化平均差, Pae et al. 2015)。

2.5〜20mg/日の用量で用量反応関係は観察されなかった(Meeker et al. 2015)。
(注: つまり効果がみられなかったので容量を増やせば効くようになるというわけではない可能性がある。もしくは少量で効果がみられる場合がある。)

最も一般的な副作用は、悪心、嘔吐、下痢および口渇である(Pae et al. 2015, Meeker et al. 2015)。

トリンテリックスは大うつ病性障害患者の認知にもプラスの効果があることが示唆されている。Rosenblat et al.(2015)によるシステマティックレビューでは、トリンテリックスが精神運動速度、実行機能、および遅延想起(SMD 0.2~0.3)の領域を改善する可能性が示されている。

これまでのところ、トリンテリックスと他の抗うつ薬との直接的な比較研究は少ない。

Llorca et al. (2014)が実施したメタ回帰分析では、トリンテリックスと、作用機序の異なる他の7つの抗うつ薬の有効性および忍容性が間接的に比較されている。

有効性に関してはこれらの薬と差異はなかった。忍容性に関しては、ベンラファキシンよりも忍容性が高く、アゴメラチンよりも忍容性が低かった。

Li et al.(2016)による系統的レビューでは、サインバルタの方がトリンテリックスよりも有効である可能性を示唆している。一方、トリンテリックスで治療した患者は副作用(例えば、悪心、下痢、多汗症、便秘、口渇、食欲減退)の発生リスクが低い。

トリンテリックスとイフェクサーを比較したRCTでは効果の差は認められなかった(Wang et al. 2015)。Montgomery et al.(2014)は、SSRISNRIで好ましい反応を得られなかった場合、トリンテリックスの方がアゴメラチンよりも有効である可能性があると述べている。

全般性不安症(GAD)

全般性不安障害へのトリンテリックスの効果は多くの研究がされており、最近出版されたPae et al. (2015)らメタアナリシスでは4つが解析の対象にされている。

  • Bidzan L, Mahableshwarkar AR, Jacobsen P, Yan M, Sheehan DV. Vortioxetine (Lu AA21004) in generalized anxiety disorder: results of an 8-week, multinational, randomized, double-blind, placebo-controlled clinical trial. Eur Neuropsychopharmacol 2012;22:847-57.
  • Mahableshwarkar AR, Jacobsen PL, Chen Y, Simon JS. A randomised, double-blind, placebo-controlled, duloxetine-referenced study of the efficacy and tolerability of vortioxetine in the acute treatment of adults with generalised anxiety disorder. Int J Clin Pract 2014a;68:49-59.
  • Mahableshwarkar AR, Jacobsen PL, Serenko M, Chen Y. A randomized, double-blind, fixed-dose study comparing the efficacy and tolerability of vortioxetine 2.5 and 10 mg in acute treatment of adults with generalized anxiety disorder. Hum Psychopharmacol 2014b;29:64e72.
  • Rothschild AJ, Mahableshwarkar AR, Jacobsen P, Yan M, Sheehan DV. Vortioxetine (Lu AA21004) 5 mg in generalized anxiety disorder: results of an 8-week randomized, double-blind, placebo-controlled clinical trial in the United States. Eur Neuropsychopharmacol 2012;22:858-66.

結果は、トリンテリックスがプラセボよりも有効であった。注目すべき薬剤は、重度のGAD患者のサンプルにおいて良好な結果(SMD=0.338)を示しているということだ。

4つのRCTすべてで最も多くみられた副作用は吐き気、頭痛、めまい、口渇であった。

容量

うつ病の患者で5〜20 mg/日。5mg/日からの投与をEU保健局は推奨している(European Medicines Agency 2013)。

半減期

トリンテリックスの血漿濃度のピークは7〜11時間(Tmax)、半減期は66時間程度。

代謝

トリンテリックスの代謝にはCYP2D6、CYP3A4 / 5、CYP2C9、CYP2C19、CYP2A6、CYP2C8、CYP2B6が関係している。CYP2D6が3-methyl-4-(2-piperazine-1-yl-phenylsulfanyl)-benzoic acid (Lu AA34443)への代謝に関連しているので、CYP2D6が主たる代謝酵素ということになる。

薬理学的特徴

このデータはSowa-Kućma et al.(2017)ではなくWikipediaから引用した。トリンテリックスの薬物のプロファイルはビラゾドンに比較的よく似ている。 ("Lundbeck's "Serotonin Modulator and Stimulator" Lu AA21004: How Novel? How Good? - GLG News". Archived from the original on 2011-07-24.)

ターゲット 親和性 機能活性 アクション
K i(nM) IC 50 / EC 50(nM) IA(%) IA (%)
SERT 1.6 5.4 阻害
NET 113 阻害
5-HT_1A 15 200 96 作動薬
5-HT_1B 33 120 55 部分作動薬
5-HT_1D 54 370 拮抗薬
5-HT_3 3.7 12 拮抗薬
5-HT_7 19 450 拮抗薬
β_1 46 拮抗薬

雑感など

以下はレビューに含まれない感想などである。あくまでも感想なので科学的な根拠は言うまでもない薄弱である。

投薬の方針

今までの抗うつ薬とは少し異なった機序が含まれるため、抗うつ薬の選択肢が広がる可能性があるのかもしれない。Montgomery et al.(2014)は、現在までSSRISNRIで十分な効果が得られなかったケースでの有効性を示唆しているので、今までの薬では改善されなかった患者への投与を検討してもよいかもしれない。

やはり、狙い目としては以前にエントリを入れた高齢者への投与であろう。
参考:トリンテリックス(ボルチオキセチン)と高齢者のうつ病

高齢者にサインバルタは有効だとはいえ、副作用で飲めないことがあったり、Poop-Out(へたばり、腰折れ)がある。長期間の投与が必要な高齢者のうつ病にはトリンテリックスはフィットする可能性がある。

追記
トリンテリックスには性機能障害の副作用が出にくいという特徴があるらしい。
この件についてエントリをした。
トリンテリックス(ボルチオキセチン)と性機能障害 - 井出草平の研究ノート  

全般性不安症

比較研究はないがSSRIの方が優れているとみて間違いないだろう。重度の全般性不安症に有効ということだが、このようなケースにはうつ病などの他の症状が併存していることが多く、全般性不安症をターゲットにしているというよりもEmotional Disorderをターゲットにしていると捉えた方がおそらく適当だろう。Emotional Disorderは昔風に言うと神経症が一番近い言葉である。

余談だがバルドキサンもSADのRCTがなく、GADのRCTがあるという似たような状況である。この種の薬剤は不安に直接作用しているとうよりも、Emotional Disorderへのアプローチをして全般的な精神健康度を上げることによって不安を低減していると考えた方が良いように思う。軽症のGADにトリンテリックスはあまり有効ではないという結果とも符合する。

半減期

血中濃度が上がるのがとても遅く、半減期がとても長い薬である。3日ほど飲み忘れても対して影響がないもしれない。

服薬量

容量依存性がないらしいので少量(5mg)でも反応する確率は十分にあるだろう。EU保健局が推奨するように5mg/日で始めるのがよいように思う。ただ、発売されるのは10mgと20mgなので、カッターで切らないと5mg錠は用意できない。EU保健局の言いつけを守るには錠剤をカッターで切断せざるを得ない。アメリカで発売される錠剤を見る限りカッターで割っても問題なさそうではある。


Vortioxetine - Wikipediaより

追記 日本で発売される錠剤の写真がウェブにアップされていた。 f:id:iDES:20191127143800j:plain

もともと半分に割る線が入っていた。アメリカで流通している錠剤より利便性が高い仕上がりとなっているようだ。

日本での薬価が高くない(10㎎錠: 168.9円/20㎎錠: 253.4円)らしい。5mgで効果が出る人は20㎎錠を4分割すると3割負担で532円/28日である。10mgでもその倍くらいである。新薬とは思えないくらい経済的である。ただ、そんなケースが多発すると、製薬会社に利益が出るだろうかと少し心配になる。

名前

アメリカではBrintellixという名前だったが、血液減量薬Brilintaと混同しやすいということで、Trintellixへ名前を変更したそうだ。EUではBrintellixのままである。(https://en.wikipedia.org/wiki/Vortioxetine#Names)

副作用の対策

副作用としてよく見られているのは悪心、嘔吐、下痢である。抗うつ薬でこれらの症状が出る場合は5-HT3レセプター関連だと相場は決まっている。SSRIなどの抗うつ薬は5-HT3作動薬なので、逆の作用を持つ拮抗薬であるイリボー(ラモセトロン)の投与をして副作用を抑制する。

今回もそのパターンかと思いきや少しイレギュラーかもしれない。SERTがメインの効果なのでSSRIと同様なので、5-HT3レセプターに作動薬として消化管症状が起こっている可能性が高いが、薬物のプロファイルをみるとトリンテリックスは、5-HT3の拮抗薬でもあり、しかもKi値は3.7と比較的強力である。どちら方向にに傾いてこれらの症状がでているか正直よくわからない。

なみに5-HT3に対してピンポイントで作動薬として働く薬はない(https://en.wikipedia.org/wiki/5-HT3_receptor#Agonists)。一覧にはないがベンゾジアゼピンには5-HT3作動薬の効果があるものの、不十分である。 5-HT3作動薬として働きやすい薬はフルボキサミンくらいだろう。副作用のために抗うつ剤を追加するのは奇妙な話なので、通常はこんなことはせずに他の薬にスイッチする。しかし、有効な薬がトリンテリックスしかなく、副作用で飲めないという場合に限っては、フルボキサミンの追加は検討してみてもよいかもしれない。

注意が必要なのは5-HT1Aレセプターも悪心、嘔吐は関係することがある点だ。フルボキサミンの追加案は5-HT3レセプター関連だと仮定しての話で、5-HT1Aが犯人であれば、フルボキサミンの追加は意味がない。

5-HT1Aの拮抗薬(https://en.wikipedia.org/wiki/5-HT1A_receptor#Antagonists)で現実的なのはβ拮抗薬である。しかし、トリンテリックスにはβ拮抗作用もあって、このあたりもよくわからない所である。

副作用としてめまいが上がっている。めまいはβ拮抗作用でよく起こる副作用である。もちろん他の作用でも起こりうるが、トリンテリックスのβ拮抗作用で発現しているならば、副作用の発現のメカニズムは少しややこしい薬だと言えよう。

ちなみに悪心・嘔吐だけであれば、プリンペランナウゼリンといったドパミン関連の薬を入れるのが通常の対応である。このパターンの副作用であればそれほど難しく考える必要はない。

ADHDへのトライアル

NETの阻害に関してのKi値が掲載されているので、ADHDに有効性があるのではと思ったところ、既に成人ADHDへの臨床研究がされていた。 https://clinicaltrials.gov/ct2/show/NCT02327013 https://www.ncbi.nlm.nih.gov/pubmed/30843450

結果はダメだったようだ。ストラテラ(アトモキセチン)のNETのKi値は5に対して、トリンテリックスは113なので、結果はやる前からわかっていたのでは?という気もしないでもない。きっとお金がどこかから降って来たのだろう。

文献

  • European Medicines Agency: BRINTELLIX (vortioxetine)1. http://www.ema.europa.eu/docs/en_GB/document_library/EPAR-Public_assessment_report/human/002717/WC500159447.pdf. In: Agency EM, editor. 2013.
  • Citrome L. Vortioxetine for major depressive disorder: a systematic review ofthe efficacy and safety profile for this newly approved antidepressant ? what is the number needed to treat, number needed to harm and likelihood to be helped or harmed. Int J Clin Pract 2014;68:60-82.
  • Baldwin DS, Chrones L, Florea I, Nielsen R, Nomikos GG, Palo W, et al. The safety and tolerability of vortioxetine: analysis of data from randomized placebocontrolled trials and open-label extension studies. J Psychopharmacol 2016;30:242-52.
  • Berhan A, Barker A. Vortioxetine in the treatment of adult patients with majordepressive disorder: a meta-analysis of randomized double-blind controlledtrials. BMC Psychiatry 2014;14:276.
  • Pae CU, Wang SM, Han C, Lee SJ, Patkar AA, Masand PS, et al. Vortioxetine: ameta-analysis of 12 short-term, randomized, placebo-controlled clinical trialsfor the treatment of major depressive disorder. J Psychiatry Neurosci2015;40:174-86.
  • Meeker AS, Herink MC, Haxby DG, Hartung DM. The safety and efficacy of vortioxetine for acute treatment of major depressive disorder: a systematic review and meta-analysis. Syst Rev 2015;4:21.
  • Fu J, Chen Y. The efficacy and safety of 5 mg/d Vortioxetine compared to placebo for major depressive disorder: a meta-analysis. Psychopharmacology (Berl) 2015;232:7-16.
  • Thase ME, Mahableshwarkar AR, Dragheim M, Loft H, Vieta E. A meta-analysis of randomized, placebo-controlled trials of vortioxetine for the treatment of major depressive disorder in adults. Eur Neuropsychopharmacol 2016;26:979-93.
  • Rosenblat JD, Kakar R, McIntyre RS. The cognitive effects of antidepressants in major depressive disorder: a systematic review and meta-analysis of randomized clinical trials. Int J Neuropsychopharmacol 2015;2015.
  • Llorca PM, Lancon C, Brignone M, Rive B, Salah S, Ereshefsky L, et al. Relative efficacy and tolerability of vortioxetine versus selected antidepressants by indirect comparisons of similar clinical studies. Curr Med Res Opin 2014;30:2589-606.
  • Li G, Wang X, Ma D. Vortioxetine versus duloxetine in the treatment of patients with major depressive disorder: a meta-analysis of randomized controlled trials. Clin Drug Investig 2016;36:509-17.
  • Wang G, Gislum M, Filippov G, Montgomery S. Comparison of vortioxetine versus venlafaxine XR in adults in Asia with major depressive disorder: a randomized, double-blind study. Curr Med Res Opin 2015;31:785-94.
  • Montgomery SA, Nielsen RZ, Poulsen LH, Haggstrom L. A randomised, doubleblind study in adults with major depressive disorder with an inadequate response to a single course of selective serotonin reuptake inhibitor or serotonin-noradrenaline reuptake inhibitor treatment switched to vortioxetine or agomelatine. Hum Psychopharmacol 2014;29:47-82.
  • Pae CU, Wang SM, Han C, Lee SJ, Patkar AA, Masand PS, et al. Vortioxetine, a multimodal antidepressant for generalized anxiety disorder: a systematic review and meta-analysis. J Psychiatr Res 2015;64:88-98.
  • Wang G, Gislum M, Filippov G, Montgomery S. Comparison of vortioxetine versus venlafaxine XR in adults in Asia with major depressive disorder: a randomized, double-blind study. Curr Med Res Opin 2015;31:785-94.

99ルール - ブートストラップは99回反復しろ!

今回は潜在クラス分析のBLRTのエントリの積み残しである。

RandomLCAパッケージを用いた潜在クラス分析とBLRT[R]

先のエントリーでは修正値を変更すればブートストラップの反復回数は連動して変えられるのでは?と書いていたが、修正値は1のままで、一般的にはブートストラップの反復回数を変化させるようである。

99ルールとはBoos(2003)で使われている言葉で99とはブートストラップの反復回数である。なぜ99回なのかというと、計算式の中で修正値として1を加算するため、99+1=100となるためである。Boosが呼んでいるだけで広く使われているわけではなさそうである。うっかり正式な場では使わないように気をつけたい。

projecteuclid.org

ブートストラップのP値はb / Bで計算する。bとは検証対象の仮説よりも、代替仮説の統計量の方が大きい数である。

潜在クラス分析でいうと、例えば2クラスモデルから3クラスモデルにしてモデルが改善したかを見る場合、3クラスモデルが検証対象の仮説になり、2クラスモデルが代替仮説になる。対数化して差を取るので尤度比の部分である。

Bはブートストラップの試行回数である。Davison & Hinkley(1997)の148ページの式を分割して書くと以下のようになる。

p_{B} =\frac{  \sharp \ T_i^* \geq T_0}{B} \tag{1}

 p_{boot} =\frac{(p_B +1)}{(B + 1)} \tag{2}

Boos(2003)によれば、一般的には(p_B +1) / (B + 1)というP値を使うそうだ。先ほどあげたブートストラップの代表的な教科書であるDavison & Hinkley(1997)でも書かれている内容である。

Bootstrap Methods And Their Application (Cambridge Series in Statistical and Probabilistic Mathematics)

Bootstrap Methods And Their Application (Cambridge Series in Statistical and Probabilistic Mathematics)

Efron & Tibshirani(1993)では異なったP値の考え方も紹介されているらしいが、それほど興味がないので、このくらいにしておこう。

An Introduction to the Bootstrap (Chapman & Hall/CRC Monographs on Statistics & Applied Probability)

An Introduction to the Bootstrap (Chapman & Hall/CRC Monographs on Statistics & Applied Probability)

潜在クラス分析のブートストラップの試行回数

こちらは潜在クラス分析のブートストラップを提案した論文に既に書いてあった。

rss.onlinelibrary.wiley.com G. J. McLachlan, 1987, "On Bootstrapping the Likelihood Ratio Test Statistic for the Number of Components in a Normal Mixture." Applied Statistics, Volume36, Issue3: 318-324.

McLachlan (1987)はシミュレーションをしてブートストラップの反復回数は最低99回必要だと述べている。

ということで、99を選択するのがよさそうである。

考察

Beathは999回を指定していたが、時間に余裕があるか、PCスペックに余裕のある時だけでよさそうだ。経験的には、反復回数を100から1000にしたからと言って、推定結果が大きく異なったということは、あまりないので、99ルールは妥当なところだと思う。

99ルールとは言っているものの、分母が割り切れる数であればいいので199回(+修正値1=200)でも問題はない。とりあえず、分数の分母に来た時に、1) 99回以上であり、2) 1を足して分母になった時に必ず割り切れる数となるようになればいいのだ。

フリーでアクセスできるものとしては、こちらの論文も参考になると思う。 europepmc.org

おまけ

先のエントリではCPUにRyzen 5 3600を搭載したPCで273.61秒と4分半くらいであった。
Ryzen 5 3600は2019年7月発売のCPUで研究者が一般的に使っているPCよりかなり処理速度が速いので、普段使っているノートPCでも同様のテストをしてみた。 CPUはIntelCore i5 7200Uである。

ベンチマークの結果は言うまでもなく、圧倒的にRyzen 5 3600が速い。 Please click the green button to continue. - UserBenchmark

f:id:iDES:20191126043835p:plain

Ryzenとは違って並列処理がされ始めた。びっくりした。

   ユーザ   システム       経過  
    396.83       1.21     424.18 

だいたい7分くらいかかっている。 性能差ほど時間差は出ていない。ユーザー時間がRyzenのシングルコア処理で272.84 だったので、並列処理がされていることも確認できる。

IntelのCPUの方がAMのCPUよりもRのブートストラップはうまく作動するのかもしれない。

少し世代は古いがSkyLake世代のPCもあるので、そちらでも計測してみた。
Core i7 6700である。Core i5 7200Uよりは古いが、性能的にはCore i7 6700の方が高い。
https://cpu.userbenchmark.com/Compare/Intel-Core-i7-6700-vs-AMD-Ryzen-5-3600/3515vs4040

f:id:iDES:20191126045848p:plain

いくつかのコアは使われているが、並列処理をしている感じではない。

   ユーザ   システム       経過  
    286.78       1.48     395.47 

6分半くらいである。ベンチマークの結果とも乖離していない結果である。こちらはユーザー時間からシングルコア処理であったことがわかる。
SkyLake世代は単純なシングルコア処理だが、KabyLake世代だとCPUの方でうまく並列処理をするようになっているのかもしれない。

ともあれRyzen 5 3600が最も速いことには違いがないので、できるだけ速いCPUでやった方がいいという、当たり前の結果ではあった。
個人的にはリモートで Ryzenを載せたPCに処理をさせようと思う。
もちろん速いCPUを使えるに越したことはないが、ここ2~3年くらいの間に作られたCoreシリーズかRyzenシリーズのPCで計算すれば問題なく使えるレベルだろう。

並列処理をするとBLRTのブートストラップは速くなるか

先日、潜在クラス分析のBLRT、要するにブートストラップを走らせるとやたらと時間がかかった。999回に設定してあることが大きな原因になっているが、Rががシングルコア処理をしているのも原因ではないかと考え、並列処理について検討してみた。MplusであればCPUのコア数を指定するコマンドがあるので、Rでもできないかと考え、試してみた。

ちなみに、先のエントリで実行したのは下記のスクリプトである。

library(randomLCA)
data(dentistry)
dentistry.lca2 <- randomLCA(dentistry[,1:5],freq=dentistry$freq,nclass=2)
dentistry.lca3 <- randomLCA(dentistry[,1:5],freq=dentistry$freq,nclass=3)
obslrt <- 2*(logLik(dentistry.lca3)-logLik(dentistry.lca2))

## 以下BLRT
nsims <- 999
thesims <- simulate(dentistry.lca2, nsims)
simlrt <- rep(NA,nsims)
for (isim in 1:nsims) {
submodel <- refit(dentistry.lca2,newpatterns=thesims[[isim]])
fullmodel <- refit(dentistry.lca3,newpatterns=thesims[[isim]])
simlrt[isim] <- 2*(logLik(fullmodel)-logLik(submodel))
}

以下BLRTの所が時間がかかる場所。

処理時間の計測

処理時間は次のようにすると計測できるようだ。標準機能である。
stat.ethz.ch

ptm<-proc.time()

スクリプト

proc.time()-ptm

処理時間を知りたいスクリプトをサンドウィッチすると計測できる。

シングルコアでの計測

上記のBLRTの部分を走らせた結果が以下。。

ユーザ   システム       経過  
 272.84       0.75     273.61

ユーザーがRが実際に動いていたCPU時間で、経過が実際の経過時間なので273.61秒ということになる。

Microsoft R Open

mran.microsoft.com

Microsoft R OpenはCPUコア数を自動で判断して並列処理をする機能がある。

Microsoft R Openについてはこちらを参照のこと。

techtarget.itmedia.co.jp

RStudioでMicrosoft R Open 3.5.3を動かすと以下のようにCPU使用率が変わった。

f:id:iDES:20191122204603p:plain

いい感じで並列処理がされている。

処理時間は以下のようになった。

ユーザ   システム       経過  
1303.42       0.62     340.08

速くなっているのかなと思ったら、遅くなっていた。あれれ?

ユーザの時間の増加はマルチコア処理がされていてCPUの使用時間の総計が多くなっているためだろう。ユーザー時間が増えて、経過時間が遅くなっているということは、並列処理をしても早くならない、計算を仕分けしたり統合したりしていると、余計に時間がかかる、ことだろうか。

doSNOW

Rでは並列化させてた処理をする方法はいくつかある。代表的なdoSNOWパッケージを利用してみた。
下記のシミュレーションでもdoSNOWは良い値を出している。

itcweb.cc.affrc.go.jp

インストールして呼び出し、処理するスクリプトを下記のようにサンドイッチする。

atm<-proc.time()
library(doSNOW)
cl.SOCK <- makeCluster(12, type="SOCK")
registerDoSNOW(cl.SOCK)

スクリプト

stopCluster(cl.SOCK)
proc.time()-atm

結果は下記。

 ユーザ   システム       経過  
  290.06       1.09     295.94

結果は悪くないが、並列化をしないよりもやや遅い。

考察

結果失敗である。
スクリプトが並列化に対応すれば、処理も早くなるのだろうが、今の形では速くはならないようだ。
ベイズ統計学をするときにも並列化は必要になるが、RStanには並列化のコマンドがデフォルトで存在しているのでMicrosoft R OpenやdoSNOWは不要である。

ブートストラップは他の分析でも使用するが、今のところ僕は潜在クラス分析の時くらいしか使わないし、人生でそれほど数を走らせるわけではないので、時間とコストを考えると何も考えずシングルコアに頑張ってもらうのでもいいのかもしれない。

poLCAでデータ加工をしてRandomLCAでBLRTの計算をする[R]

私たちが使い慣れているデータ(行にケース、列に変数という形式のもの)を使ってRでブートストラップ尤度比検定(Bootstrapped Likelihood Ratio Test)を計算する方法についてのエントリである。
poLCAには残念ながらBLRTの機能がないので、poLCAのpredcell関数でデータを加工してからRandomLCAパッケージでBLRTの計算をするという手順を踏む。

データ

poLCAのcarcinomaデータを利用する。 子宮頸部のがんのデータである。デモなのでデータの説明は省く。

library(poLCA)
data(carcinoma)
carcinoma[100:105,] #100~105行目を抽出

次のようなデータである。

    A B C D E F G
100 2 2 2 2 2 1 2
101 2 2 2 2 2 1 2
102 2 2 2 2 2 1 2
103 2 2 2 2 2 2 2
104 2 2 2 2 2 2 2
105 2 2 2 2 2 2 2

predcell関数でデータ形式を変換

データ形式を下記の手順で変換する。

data(carcinoma)
f <- cbind(A,B,C,D,E,F,G)~1
clc2 <- poLCA(f,carcinoma,nclass=2)
carcinoma2<- clc2$predcell # スタック形式carcinoma2とする
print(carcinoma2) 

このようにするとRandomLCAで扱える形式のデータに変換できる。

   A B C D F G observed expected
1  1 1 1 1 1 1       36   30.028
2  1 2 1 1 1 1       10   16.197
3  1 2 1 1 1 2        6    2.023
4  2 1 1 1 1 1        2    3.751
5  2 1 2 1 1 2        1    0.203
6  2 2 1 1 1 1        4    2.023
7  2 2 1 1 1 2        8    4.075
8  2 2 1 1 2 2        1    2.769
9  2 2 1 2 1 2        3    4.447
10 2 2 1 2 2 2        3    3.222
11 2 2 2 1 1 2       13   11.858
12 2 2 2 1 2 2        5    8.592
13 2 2 2 2 1 2       10   13.797
14 2 2 2 2 2 2       16    9.997

もっとスマートな方法は他にもありそうだが、とりあえずpredcell関数しか僕は発見できていないので、今回はこの方法を利用する。

randomLCAパッケージで潜在クラス分析

データcarcinoma2のA~Gはデータは1/2で構成されている。RandomLCAは0/1である必要があるため、memiscパッケージを利用して0/1にリコードする。

library(memisc)
carcinoma2$A <- recode(carcinoma2$A , 0 <- 1, 1 <- 2)
carcinoma2$B <- recode(carcinoma2$B , 0 <- 1, 1 <- 2) 
carcinoma2$C <- recode(carcinoma2$C , 0 <- 1, 1 <- 2) 
carcinoma2$D <- recode(carcinoma2$D , 0 <- 1, 1 <- 2) 
carcinoma2$E <- recode(carcinoma2$E , 0 <- 1, 1 <- 2) 
carcinoma2$F <- recode(carcinoma2$F , 0 <- 1, 1 <- 2) 
carcinoma2$G <- recode(carcinoma2$G , 0 <- 1, 1 <- 2) 

データタイプは整数型でも実数型でも走るようだ。

今回は2クラスモデルと3クラスモデルを比較する。
比較する2クラスモデルと3クラスモデルを作成しておく。

library(randomLCA)
carcinoma2.lca2 <- randomLCA(carcinoma2[,1:7],freq=carcinoma2$observed,nclass=2)
carcinoma2.lca3 <- randomLCA(carcinoma2[,1:7],freq=carcinoma2$observed,nclass=3)

BLRT

BLRTのスクリプトに関する細かい説明は先のエントリを参照のこと。
まず、対数尤度の差の二乗。

## 3クラス、2クラスモデルの対数尤度の差の二乗
obslrt <- 2*(logLik(carcinoma2.lca3)-logLik(carcinoma2.lca2)) 
print(obslrt)

結果は以下。

'log Lik.' 47.09539 (df=23)

BLRT。

nsims <- 999 # ブートストラップの回数を指定
thesims <- simulate(carcinoma2.lca2, nsims) # シミュレーションを作成
simlrt <- rep(NA,nsims) 
for (isim in 1:nsims) {
submodel <- refit(carcinoma2.lca2,newpatterns=thesims[[isim]]) # submodelは比較対象の2クラスモデル
fullmodel <- refit(carcinoma2.lca3,newpatterns=thesims[[isim]]) # fullmodelは3クラスモデル
simlrt[isim] <- 2*(logLik(fullmodel)-logLik(submodel)) # 2つのモデルの尤度の差を二乗
}

P値。

print((sum(simlrt>=obslrt)+1)/(nsims+1)) # P値の計算

結果は以下。

p-value=0.001

Mplusに比べて手順が煩雑なのが欠点であろうか。