井出草平の研究ノート

アレン・フランセス「臨床的意義の重要性」

www.psychiatrictimes.com

July 2, 2010 Allen Frances, MD

特徴的なクラスターの有無だけで精神障害を定義してはどうでしょうか?なぜ苦痛や障害も必要だと感じたのでしょうか?

精神症状は一般集団ではかなり広範囲にみられ、ほとんどの正常人には少なくとも1つ、多くの人にはいくつかみられます。単独で存在する場合、1つの症状 (または少数の症状) で精神疾患が生じることはありません。症状が精神障害の一部とみなされるには、さらに2つの条件を満たす必要があります。第1に、関連するDSM基準セットに示されているように、特徴的な症状のクラスターが存在しなければならない。うつ病、不安、不眠、記憶障害、注意力の問題などの単独の症状だけでは、精神障害の診断を正当化するのに十分ではありません。

第二、ここでの主な対象は、症状が臨床的に有意な苦痛または社会的または職業的機能における臨床的に有意な障害を引き起こすことです。この注意事項は非常に重要であるため、ほとんどのDSM基準セットには別の項目として含まれています。症状があるだけでは十分ではありません。生活に深刻な問題を起こしていなければなりません。

特徴的なクラスターの有無だけで精神障害を定義してはどうでしょうか?なぜ苦痛や障害も必要だと感じたのでしょうか。DSM障害のほとんどは、重症度の段階的スペクトルに沿って存在する。重症になると、症状から生じる患者の苦痛や障害が非常に明白になり、その症状が精神障害であることを疑う余地はありません。しかし、ほとんどの障害の軽症のものでは、正常と精神障害を区別する明確な境界線はありません。

何が臨床的に重要なのかをどのように定義すればよいのでしょうか。残念ながら、これは正確な指標を持たない曖昧な言葉なのです。臨床的に重要な精神障害を有するのに十分な苦痛または障害を経験しているかどうかを判断することは、本質的に困難で主観的な判断であり、客観的な基準なしに行われなければなりません。

役立つヒントをいくつか紹介します。まず、明確な正解はないことを理解しましょう。自分の状態が診断と治療を正当化するほど重度であるかどうかの質問に答えるには、少なくともいくらかの不確実性が不可避であることを受け入れなければならない。この認識は、いくつかの重要な意味をもたらします。注意深く待つことは、どちらか一方の結論に飛びつくよりも、最初の一歩としてはずっと良いかもしれません。特に、それほど長くは続かず、障害もないような軽い症状には、副作用のない安上がりな治療薬がよく使われます。

次に、これは主観的すぎて自分では決められないかもしれません。ある人はストイックで、自分が困っていることを受け入れる前に、文字通り死の淵に立たされなければならないでしょう。一方、完璧主義者、感情的な心気症患者、破局論者は、日常生活で予想される痛みに過ぎないかもしれないのに、診断と治療を求めたりします。疑わしいと感じたら、セカンドオピニオン、サードオピニオン、フォースオピニオンを得るようにしましょう、臨床家からも、知識のある家族からも。

このような状況では、通常、臨床医の助けを借りて、診断を受けることのプラス面とマイナス面のリスク/ベネフィット評価を行うことが有効です。基本的な質問は「この診断を受けることは、私を助ける可能性が高いか、それとも私を傷つける可能性が高いか」に尽きます。安全で効果的であることが証明されている治療法がある場合は、その診断を受け入れることが賢明ですが、証明された治療法がない場合、または利用できる治療法に危険な副作用や合併症がある場合は、疑わしい診断は避けなければなりません。

臨床医は、診断の有無が明確でない場合でも、治療の試行を提案することがある。その理由は、治療によって気分が良くなるのであれば、診断基準を完全に満たしているかどうかを気にする必要はない、あるいは、肯定的な反応が得られれば、その治療がおそらく必要であることが証明される、というものです。これらはいずれも不正確で誤解を招く議論です。軽度の疾患では、プラセボの反応率が非常に高く、しばしば約50%であり、これは薬物療法で得られる反応率に非常に近いものです。軽度の疾患で薬を飲み始めると、何が原因で症状が改善したのか、プラシーボ効果なのか、薬の有効成分なのか、わからなくなります。また、薬の効果で改善したと勘違いして、むやみに薬を長く飲み続けることになり、副作用のリスクもあります。ですから、軽い症状の場合は、まず時間をかけて注意深く待つこと(このパッケージには常に運動療法が含まれます)、次に心理療法、そして薬物療法は3番目と最後の手段という順番がベストでしょう。特に、その診断が、抗精神病薬の使用など、かなり危険な治療の引き金になる可能性がある場合には、プライマリーケア医に5分ほど診てもらっただけで薬をもらうのは、通常良い考えとは言えません。

あなたの精神状態が「中間的」である場合、つまり、すべてが順調というわけではないが、明らかに障害があるというわけでもない場合、診断や治療に関してすぐに決断してはいけないということです。しばらく時間をおいて、他の人の意見も聞いてみましょう。通常、1〜2ヶ月で事態は収拾します。症状が改善されれば、診断の必要はありませんでしたが、症状が長引いたり、悪化したりした場合は、助けを求めてください。

傾向スコアによる重み付けを行うPSweightパッケージ[R]その2

こちらの続き。

ides.hatenablog.com

ides.hatenablog.com

cran.r-project.org

こちらの4章の翻訳。

https://cran.r-project.org/web/packages/PSweight/vignettes/vignette.pdf


4.3. 複数の処置法を用いた傾向スコアによる重み付け

セクション4.2で二重処置のケースで提供した構文は,多重処置のケースに継ぎ目なく一般化することができる。したがって,冗長にならないように,このサブセクションの目的は,同じ手順を繰り返さないことである。その代わりに,複数の治療を伴うPSweightの追加的な特徴を指摘することによって,最後のサブセクションを補完する。シンプルにするために,IPWと重複を改善する3つのタイプの重みに焦点を当てることにする. OW、MW、EW(Li and Li 2019)である。

一般化傾向スコアの推計とバランス評価

3水準変数であるDmultを治療対象として使用する。 人口の約半数が高度な学歴を有し、中程度の学歴を有する者と学歴を有しない者がほぼ同数である。また、比率推定量の推定と推論を説明するために、賃金の二値化結果であるwagebinを用いる。この二値化された賃金は、19913年の30-39歳のイギリス人男性労働者の平均時給のカットオフ値として得られた。平均時給は8.23円であり、log(8.23) ≒ 2.10をカットオフとした。平均時給を上回る者は1610人、下回る者は2032人であり、平均時給を上回るための学歴の一対(加重)平均治療効果を推定することに関心がある。
一般化傾向スコアを推定するために,二値治療の場合と同じ共変量セットを用いて,多項式回帰モデルps.multを指定する.

ps.mult <-  Dmult  ~  white + maemp + as.factor(scht) + as.factor(qmab) +
  as.factor(qmab2) + as.factor(qvab) + as.factor(qvab2) + paed_u   + maed_u +
  agepa   + agema + sib_u + paed_u*agepa   + maed_u*agema

そして、傾向スコア推定値を得て、SumStat()関数で重み付き共変量バランスを評価する。 このコンポーネントは、一般化傾向スコアを可視化するために密度プロットのみを許可する(ヒストグラムは許可しない)ことを除いて、2値処置のケースと同様である.具体的には、type = "hist "を指定しても、plot.SumStat()関数は密度プロットを返す。この場合,"Histogram only available for binary treatment "という警告メッセージが生成される。

bal.mult <-  SumStat(ps.formula = ps.mult,
    weight = c("IPW", "overlap",  "matching", "entropy"),
    data = NCDS) 
plot(bal.mult,  type = "hist")

Warning  message:
Histogram only available for binary treatment. Density plot provided instead.

一般化傾向スコアの分布を図3に示す(治療群名のアルファベット順)。 advanced qualification (">=A/eq") または no qualification ("None") を受ける一般化傾向スコアについては、グループ固有の分布が分離しているため、重複が軽度に欠如していることがわかる。bal.multには4つの重み付けスキームが含まれているので、最大ペアワイズASDをプロットし、(重み付け)共変量バランスを1つのLoveプロットで評価かる。

plot(bal.mult,  metric = "ASD")

共変量は重み付けをする前に3つのグループ間でアンバランスになっている。IPWは一般的に共変量バランスを改善することができまる、重複の不足により最大ペアワイズASDはまだ時々閾値0.1を超えることがある。一方、OW、MW、EWはいずれも重複が改善された部分母集団を強調し、すべての共変量にわたってより良いバランスを提供する。

一般化傾向スコアのトリミング

PSweightパッケージは、(一般化された)傾向スコアに基づいてトリミングを行うことができる。IPWは、図4の3つのグループ間の共変量を十分にバランスさせないので、IPWのバランスを改善する方法としてトリミングを検討する。PSweightパッケージによって実行されるトリミングには、(1)極端な(一般化傾向スコア)を持つユニットを削除する対称トリミング(Crump et al.2009;Yoshida et al. 2018)と、(2)(ペアワイズ)ATEを推定するための最もffiient IPW定量を提供する最適トリミング(Crump et al.2009; Yang et al. 2016)の2種類がある。具体的には、対称的なトリミングは、デルタ引数を通じてSumStat()とPSweight()関数の両方によってサポートされている。両関数とも、Li et al. (2019)の推奨に従ってトリミング後に(一般化)傾向スコアモデルを再フィットする。また、対称トリミングと最適トリミングの両方を実行するスタンドアロンPStrim関数も提供する。 Yoshida et al. (2018)に従い、3つの治療群で、推定された一般化傾向スコアがδ = 0.067より小さいすべての個体を除外する。この閾値は、advanced qualification グループのかなりの量の個体を除去する(情報は、SumStatオブジェクトのtrim要素から引き出すことができる)。 Yoshida et al. (2018)で議論されているように、傾向トリミングはATEとATTの推定を改善することができるが、ATOとATMの推定にはほとんど影響を及ぼさない。明らかに、図5は、トリミングされたサンプルにおいて、IPWがすべてのペアワイズASDを10%以内に制御していることを示している。また、OW、MW、EWについては、トリミングによる加重バランスへの影響はほとんどなかった。

bal.mult.trim <-  SumStat(ps.formula = ps.mult,
        weight = c("IPW", "overlap",  "matching", "entropy"),
        data = NCDS, delta = 0.067)
bal.mult.trim
1050 cases trimmed,  2592 cases remained 

trimmed result by trt group: 
         >=A/eq None O/eq
trimmed     778   71  201
remained   1028  824  740

weights estimated for:  IPW overlap matching entropy 
plot(bal.mult.trim,metric = "ASD")

図5:δ=0.067の対称トリミング後の、最大ペアワイズASDメトリックを用いた3レベル処理変数Dmultのラブプロット。このプロットは、PSweightパッケージのplot.SumStat()関数によって生成されたもの。

また、トリミングの閾値を指定しない場合、PStrim関数はデータに基づいて最適な閾値を特定する最適トリミング手順をサポートする。 構文の例を以下に示す。トリミングの要約統計量を引き出すと、最適トリミングでは、上級資格者、中級資格者、無資格者において、それぞれ27%、9%、2%が除外されていることがわかる。 この排除はδ=0.067の対称的トリミングと比較してより保守的である。しかし、最適トリミング後の共変量バランスは、図5と同様であり、省略した。

PStrim(ps.formula = ps.mult, data = NCDS, optimal = TRUE)
Summary of the trimming:
         >=A/eq None O/eq
trimmed     479   21   82
remained   1327  874  859

一組の(加重)平均治療効果の推定と推論

2.5節で紹介した比率推定量をバイナリーアウトカムのwagebinを用いて推定する。ここでは、説明のために、トリミングを行わないデータに基づいて因果効果を推定するのみであり、トリミングしたデータによる解析も全く同じ手順で行う。 多項ロジスティック傾向スコアモデルに基づいて、IPWにより結合母集団における一対の因果関係RRを求める。

ate.mult <-  PSweight(ps.formula = ps.mult, yname = "wagebin", data  = NCDS,
               weight = "IPW")
contrasts.mult <-  rbind(c(1,-1, 0), c(1, 0,-1), c(0, -1, 1))
sum.ate.mult.rr <-  summary(ate.mult, type = "RR",  contrast = contrasts.mult) 
sum.ate.mult.rr
Closed-form inference: 

Inference in log scale: 
Original group value:  >=A/eq, None, O/eq 

Contrast: 
           >=A/eq None O/eq
Contrast 1      1   -1    0
Contrast 2      1    0   -1
Contrast 3      0   -1    1

            Estimate Std.Error       lwr     upr  Pr(>|z|)    
Contrast 1  0.607027  0.115771  0.380120 0.83393 1.577e-07 ***
Contrast 2  0.459261  0.052294  0.356767 0.56176 < 2.2e-16 ***
Contrast 3  0.147766  0.121692 -0.090746 0.38628    0.2246    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

適切な対比行列を用意し、summary.PSweight()関数で対数スケールでの平均値の比較を行い、対比ベクトルaに対して \lambda^{h, R R}(a)を推定する。p値は、弱い因果関係の帰無仮説 H_0: \lambda^{h, R R}(a) に対する統計的証拠を提供する。その結果、複合母集団において、全員が高学歴の場合、平均以上の時給を受け取る割合は、全員が無学歴の場合のexp(0.607)=1.83倍となることが分かった。また、全員が上級資格を取得したときに平均以上の時給を受け取る割合は、全員が中級資格を取得したときのexp(0.459)=1.58倍である。両効果とも0.05水準で有意であり、対応する因果の帰無に対して強い証拠を与える(p値<0.001)。しかし、全員が中級の資格を取得した場合、平均以上の時間給を受け取る割合は、資格のない場合と比較してわずかに高くなるだけであり、p値は0.05を超えた。因果関係RRとその信頼区間を直接報告するには、summary.PSweight()関数によって提供される点推定値と信頼区間を単純に指数化すればよい。

exp(sum.ate.mult.rr$estimates[,c(1,4,5)])
           Estimate       lwr      upr
Contrast 1 1.834968 1.4624601 2.302358
Contrast 2 1.582904 1.4287028 1.753749
Contrast 3 1.159241 0.9132496 1.471493

観察された共変量において最も重複する標的集団に焦点を当てて,著者らはさらに対の因果的RRを推定するためにOWを使用した。OWは、対比較のための最良の内部妥当性を理論的に提供する;図5はまた、OWが重複集団間のより良好な共変量バランスを達成することを示す。サマリーによって提供される結果の指数。PSweight () 関数を用いて,各ペアワイズ因果RRが重複加重集団の中で大きな効果サイズを持つことを観察した。興味深いことに、重複集団の中では、誰もが中間的な資格を取得したときに平均以上の時間給を受け取る割合は、誰もが学術的な質を獲得していないときの約1.55倍になり、関連する95% CIは無効を除外する。更に,一対比較に対する標準誤差は, OW対IPWを用いた場合,より小さく, OW分析は,等平衡を有するポピュレーションに焦点を絞ることにより,一般的に増加した力に相当することを示した。MWとEWの両方を使用して分析を繰り返す;この解析の結果はOWに類似しているため、簡略化のため省略している。

ato.mult <-  PSweight(ps.formula = ps.mult, yname = "wagebin", data  = NCDS,
               weight = "overlap")
sum.ato.mult.rr <-  summary(ato.mult, type = "RR",  contrast = contrasts.mult)
exp(sum.ato.mult.rr$estimates[,c(1,4,5)])
           Estimate      lwr      upr
Contrast 1 2.299609 1.947140 2.715882
Contrast 2 1.527931 1.363092 1.712705
Contrast 3 1.505048 1.257180 1.801785

上記の結果は、重複集団の中で、高度な資格認定と中間的な資格認定とを比較した因果RRは、中間的な資格認定と資格なしとを比較した場合と同様の大きさであることを示唆している。帰無仮説 H_{0}:  \mu_3^h / \mu_2^h =  \mu_2^h / \mu_1^h  (セクション2.5も参照)に基づいて、2つの連続した原因RRの等価性を形式的にテストすることができる。操作上、対応するコントラストベクトルcontrast =c(1、1、-2)を指定する必要がある。この帰無仮説ほテストするためのp値は0.91(簡潔にするために出力を省略)であり、0.05レベルでの連続的因果RRの同等性に対する証拠の欠如を示唆する。

summary(ato.mult, type = "RR", contrast = c(1, 1, -2), CI = FALSE)
Closed-form inference: 

Inference in log scale: 
Original group value:  >=A/eq, None, O/eq 

Contrast: 
           >=A/eq None O/eq
Contrast 1      1    1   -2

           Estimate Std.Error z value Pr(>|z|)
Contrast 1  0.01509   0.12821 0.11769   0.9063
sum.ato.mult.or <-  summary(ato.mult, type = "OR",  contrast = contrasts.mult) 
exp(sum.ato.mult.or$estimates[,c(1,4,5)])
           Estimate      lwr      upr
Contrast 1 3.586050 2.841383 4.525879
Contrast 2 2.050513 1.696916 2.477791
Contrast 3 1.748855 1.375483 2.223578

最後のステップとして, OWを結果回帰と組み合わせる方法を説明し,重複集団間のペアワイズ因果RRを推定した。4.2節と同様に、二項結果回帰モデルでも同じ共変量の集合を用いる。

out.wagebin <-  wagebin   ~  white + maemp + as.factor(scht) + as.factor(qmab) +
           as.factor(qmab2) + as.factor(qvab) + as.factor(qvab2) + paed_u   + maed_u +
           agepa   + agema + sib_u + paed_u * agepa + maed_u *agema

この結果の回帰式をPSweight () 関数にロードし、family="binomial"を指定して結果の種類を指定すると、logRRスケールで拡張された重複重み付け推定値が得られる。点推定値と信頼限界を乗じると、1対の因果的RRが報告される。拡張OW推定器により報告されたペアワイズ因果RRは単純OW推定器により報告されたものと類似している;さらに、信頼区間の幅もアウトカム増強の前後で同等であり、ペアワイズRRに基づく因果関係の結論は変わらない。単純なOW推定器と拡張OW推定器の間の類似性は, OWそれ自身が既に効率的であるかもしれないことを意味する。

ato.mult.aug <- PSweight(ps.formula = ps.mult, yname = "wagebin", data = NCDS,
          augmentation = TRUE, out.formula = out.wagebin, family = "binomial")
sum.ato.mult.aug.rr <- summary(ato.mult.aug, type = "RR",
          contrast = contrasts.mult)
exp(sum.ato.mult.aug.rr$estimates[,c(1,4,5)])
           Estimate      lwr      upr
Contrast 1 2.310628 1.957754 2.727105
Contrast 2 1.540176 1.375066 1.725111
Contrast 3 1.500237 1.253646 1.795331

4.4. 機械学習を使用した傾向スコアと潜在的な成果の推定

既定の一般化された線形モデルの代わりに、より高度な機械学習モデルを使用して傾向スコアと潜在的な結果を推定できる。柔軟な傾向スコアと結果推定は,モデルのミス・スペシフィケーションによるバイアスを低減し,共変量バランスを潜在的に改善することが示されている (Lee, Lessler and Stuart 2010;Hill 2011;McCaffrey et al.2013) 。これは、この方法を一般化ブーストモデル (GBM)またはSuperLearnerの方法として指定することにより、バランスチェックと構造化重み付き推定器の両方に対してPSweightで実現できる。これらのメソッドの追加のモデル仕様は、ps.controlとout.controlで提供できる。gbmにもSuperLearnerにも含まれていない機械学習モデルは、外部から推定し、ps.estimateおよびout.estimate引数を使用してインポートできる。これらの2つの議論は、傾向スコアの外部からの推定値とpotential outcomeモデルを容易に組み込むことができる場合に、PSweightの有用性を広げる。

ここでは、デフォルトの一般化線形モデルの代替としてGBMを使用する方法について説明する。この図は、二項教育 「Dany」に基づいている。GBMは,予測因子と転帰(Friedman、Hastie、Tibshirani 2000)の間の柔軟な非線形関係を可能にするノンパラメトリックツリーに基づく回帰のファミリーである。次の傾向モデル式が指定されている。ブースト回帰はすでに非線形効果と相互作用(McCaffrey, Ridgeway, and Morral (McCaffrey, Ridgeway, and Morral 2004))を捉えることができるため、式には相互作用項は含まれない。この図では、AdaBoost (Freund and Schapire 1997) アルゴリズムを使用して、制御設定ps.control=list (distribution="adaboost") によって傾向モデルを適合させている。ツリー数 (n.trees=100) 、相互作用深さ (interaction.depth=1) 、終端ノードの最小観測数 (n.minobsinnode=1) 、収縮低減 (shrinkage reduction=0.1) 、および袋詰め率 (shrinkage=0.5)など、他のモデルパラメータには既定値を使用する。これらのパラメーターの代替値は、ps.controlを介して渡すこともできる。

ps.any.gbm <-  Dany ~  white + maemp + as.factor(scht) + as.factor(qmab) +
           as.factor(qmab2) + as.factor(qvab) + as.factor(qvab2) + paed_u  +
           maed_u+ agepa   + agema + sib_u
bal.any.gbm <-SumStat(ps.formula  = ps.any.gbm, data= NCDS, weight = "overlap",
          method  = "gbm", ps.control = list(distribution = "adaboost"))

plot.SumStat()`によるバランスチェックでは、重み付け後のすべての共変量のSMDが0.1以下となり、共変量バランスが大幅に改善されたことが示された。 バランスを評価し、傾向スコアモデルの妥当性を確認した後、我々はさらにデフォルトのロジスティック回帰とパラメータを用いて、GBMを用いてアウトカムモデルを推定する。 PSweight()関数では、ps.method = "gbm" と out.method = "gbm "の両方を指定し、out.control引数はデフォルトのままにしておくことができる。 詳細なコードと出力の要約は以下の通りである。 GBMはデフォルトで共変量間の相互作用を考慮するため、ここでは相互作用項を用いずに傾向スコアモデルを再定義している。 GBMを用いた結果は、一般化線形モデルを用いた結果と非常によく似ている。

out.wage.gbm <- wage ~ white + maemp + as.factor(scht) + as.factor(qmab) +
  as.factor(qmab2) + as.factor(qvab) + as.factor(qvab2) + paed_u +
  maed_u + agepa + agema + sib_u
ato.any.aug.gbm <- PSweight(ps.formula = ps.any.gbm, yname = "wagebin",
  data = NCDS, augmentation = TRUE, out.formula = out.wage.gbm,
  ps.method = "gbm", ps.control = list(distribution = "adaboost"),
  out.method = "gbm")
summary(ato.any.aug.gbm, CI = FALSE)
Using 100 trees...

Using 100 trees...


Closed-form inference: 

Original group value:  0, 1 

Contrast: 
            0 1
Contrast 1 -1 1

           Estimate Std.Error z value  Pr(>|z|)    
Contrast 1 0.187532  0.018462  10.158 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

5. 概要

傾向スコア重み付けは因果推論と比較有効性研究のための重要なツールである。本論文では、PSweightパッケージを紹介し,二値および複数の治療群の状況でNCDSデータ例を用いてその機能を示した。PSweightは、観察研究のデザインを支援するための読みやすいバランス統計とプロットを提供するだけでなく、加法スケールと比率スケールの両方に対する (加重された) 平均処理効果について、様々な加重スキームを用いた点と分散の推定も提供する。これらの重み付けスキームはLiら (2018) 及びLiとLi (2019) で最近導入された最適重複重みを含み,等姿勢での集団間の有効な因果比較有効性証拠の生成を助けることができた。

PSweightパッケージは、傾向スコア重み付け分析のための他の有用なコンポーネントを含めるために継続的に開発されている。具体的には、PSweightの将来のバージョンには、バランシングウェイトと柔軟な変数選択ツール(Yang, Lorenzi, Papadogeorgou, Wojdyla, Li, and Thomas 2020)を使用して、事前に指定されたサブグループ分析を可能にするコンポーネントが含まれる。また、イベントまでの時間の結果と複雑な調査デザインを持つ重複重み付け推定器を研究している。これらの新機能は,方法論の拡張と同時に活発に開発されている。

傾向スコアによる重み付けを行うPSweightパッケージ[R]その1

cran.r-project.org

こちらの4章の翻訳。

https://cran.r-project.org/web/packages/PSweight/vignettes/vignette.pdf


4. NCDSデータによるケーススタディ

NCDS(National Child Development Survey)データを用いて、教育達成度の時間給への因果関係を推定するケーススタディにおいて、PSweightを実証する。4.1 節では、研究の概要を説明する。4.2節と4.3節では、それぞれ二値処理と三値処理の文脈における平均処理効果の傾向スコアによる重み付け分析を行う。セクション4.4では、傾向スコアと潜在的な結果に対する機械学習法がPSweightにどのように実装されるかを示す。

4.1. NCDS データ概要

NCDS(The National Child Development Survey)は、1958年にイギリス(UK)で生まれた子どもを対象とした縦断的な研究である(https://cls.ucl.ac.uk/cls-studies/1958-national-child-development-study/)。NCDSは17,415人の教育達成度、家族的背景、社会経済的・健康状態などの情報を収集した。Battistin and Sianesi (2011)に従って前処理を行い、1991 年に雇用された男性 3,642 人のうち、学歴と賃金情報が完全なサブセットを得て、分析を行った。 また、Multiple Imputation by Chained Equations (MICE, Buuren and Groothuis-Oudshoorn 2010) を用いて、欠損共変数をインプットし、1つのインプットデータセットを得て、その後のすべての分析に用いている*。

*我々が検討した治療前の共変量12個のうち10個に欠損Missingnessがある。最も小さい欠損の割合は4.9%で、最も大きい欠損の割合は17.2%である。我々は、説明のために1つの代入された完全なデータセットを考慮したが、より厳密な分析は、Rubinのルールによって複数の代入れたデータセットからの分析を組み合わせることによって進められることに留意されたい。

治療変数は教育達成度である。処理変数は学歴であり、二項治療の場合、学歴を取得したかどうかを示すDanyを作成した。学歴がある者は2399人であり、学歴がない者は1,243人である。複数回治療の場合、Dmultは3つのレベルで作成した。">=A/eq"と "O/eq"と "None"である。上級資格者(1,806人)、中級資格者(941人)、無資格者(895人)を表している。前処理の12の共変量または潜在的な交絡因子を考慮した。変数whiteは白人種であると自認しているかどうか、schtは16歳時に通っていた学校の種類、qmabqmab2は7歳と11歳の時の数学テストの得点、qvabqvab2は7歳と11歳の時の読解テストの点数、sib_uは兄弟の数、agepaagemaは1974年の両親の年齢、同年には母親maempの雇用形態も収集、paid_umaid_uは両親の教育年数、agepaagemaは1974年の親の年齢である。

研究変数に関する情報は、str()関数を用いて以下のように要約することができる。

library(PSweight)
str(NCDS)
'data.frame':    3642 obs. of  16 variables:
 $ white  : int  1 1 1 1 1 1 1 1 1 1 ...
 $ wage   : num  2.57 2.04 1.72 2.2 2.48 ...
 $ Dany   : int  1 1 0 1 1 0 0 1 1 1 ...
 $ Dmult  : chr  ">=A/eq" ">=A/eq" "None" "O/eq" ...
 $ maemp  : int  0 0 0 0 1 1 0 1 0 1 ...
 $ scht   : int  2 1 1 3 1 2 1 1 1 3 ...
 $ qmab   : int  2 5 4 5 3 1 4 5 5 2 ...
 $ qmab2  : int  2 5 4 4 3 1 1 4 4 5 ...
 $ qvab   : int  1 5 4 4 2 2 2 4 3 4 ...
 $ qvab2  : int  2 5 5 5 3 2 1 3 1 5 ...
 $ paed_u : int  9 0 0 10 9 10 0 11 10 10 ...
 $ maed_u : int  9 0 0 10 9 9 0 11 9 10 ...
 $ agepa  : int  60 56 57 40 57 43 43 46 43 47 ...
 $ agema  : int  59 56 53 41 45 42 38 45 43 40 ...
 $ sib_u  : int  3 0 0 1 1 1 1 1 0 3 ...
 $ wagebin: num  1 0 0 1 1 0 1 1 1 1 ...

4.2. 二元的処置による傾向スコア重み付け

傾向スコアの推計とバランスチェック

学歴の取得が時間当たり賃金の上昇につながるかどうかの因果関係を推定したいとする。学歴の取得は無作為化されておらず、潜在的な交絡因子によって影響を受ける可能性があるため、以下の傾向スコアモデルを設定し、重み付け分析を行う。

ps.any <- Dany ~  white + maemp + as.factor(scht) + as.factor(qmab) +
          as.factor(qmab2) + as.factor(qvab) + as.factor(qvab2) + paed_u + maed_u +
          agepa + agema + sib_u + paed_u * agepa + maed_u * agema

共変量の主効果に加えて,Battistin and Sianesi(2011)に従って,親の年齢と教育の交互作用を調整することを検討した。Sumstat()関数を用いてロジスティック傾向スコアモデルを推定し、IPW、処理重み、OWの3種類の重み付けスキームにおけるバランス統計量を求めた。

bal.any <- SumStat(ps.formula = ps.any, data = NCDS,
             weight = c("IPW", "overlap", "treated"))

Sumstat() 関数からの画面上の出力は、重みの選択と、重み引数に "treated" が含まれる場合にのみ選択される治療グループ (trtgrp) である。 この例では、trtgrpが指定されていないので、Sumstats()は自動的に処理変数のアルファベット順の最後のレベルを処理グループとする。Dany = 1.

bal.any
trt group for PS model is:  1 
weights estimated for:  IPW overlap treated

SumStatは,治療群レベル(ATTを定義するため)("trtgrp"),推定傾向スコア("propensity"),各重み付けスキームにおける推定重み("ps.weights"),有効サンプルサイズ("ess"),各重み付けスキームにおけるバランス統計(例: "unweighted.sumstat", "IPW.sumstat", "overrap.sumstat", "treated.sumstat")などを含むリストとして返される。さらに、各重み付けスキームのバランス統計は、共変量の非重み付け標準偏差または重み付け標準偏差の両方と、ASDとPSDの両方を含んでいる。 plot.SumStat()関数は、推定傾向スコアと共変量balance統計量の分布を可視化する。引数type = "hist"を指定すると、処置(trtgrpで定義される治療)を受ける推定傾向スコアのヒストグラムが生成される。

また、type = "density"は、各治療水準を受ける推定確率の密度を示す。図1は、我々の分析における推定傾向スコアのヒストグラムと密度プロットである。ヒストグラムは、2つのグループの軽微な分離のため、重複がわずかに欠けている可能性を示唆している。

plot(bal.any,  type = "density")
plot(bal.any,  type = "hist")

図1:plot.SumStat()関数で生成した二元処置変数Danyに関する推定傾向スコアのヒストグラムと密度プロット

最後に、plot.SumStat()で引数type = "balance"を指定すると、ASDメトリック(metric = "ASD")または最大PSDメトリック(metric = "PSD")に基づいたラヴプロットが生成されます。図2は、重み付き標準偏差(デフォルトではweighted.var = TRUE)を用いたPSDベースのラブプロットを示している。明らかに、重み付けされていない平均の乖離は、一般的に使用されるバランス閾値0.1よりも大幅に大きく、一方で傾向スコア重み付けは一般的に共変量バランスを改善することができる。 3つの重み付けのうち、OWとIPWは各共変量の最大PSDが0.1以下になるように制御されており、OWは各共変量の最大PSDが0に近く、最もバランスが取れていると言える。

plot(bal.any, type = "balance", metric = "PSD")

図2:PSweightパッケージのplot.SumStat()関数で生成した、最大PSDメトリックを用いたバイナリ処理変数Danyのラブプロット

加重平均治療効果の推定と推論

図2によれば、IPW、治療重み、OWは適切なバランスを保っているので、これら3つの重み付けスキームを用いて、ATE、ATT、ATOを推定する。傾向スコアモデルps.anyに基づき、まずPSweight()関数を用いてIPWを用いた結合集団の平均潜在的転帰を求める。

ate.any <-  PSweight(ps.formula = ps.any, yname = "wage", data  = NCDS,
                    weight=  "IPW")
ate.any
Original group value:  0, 1 

Point estimate: 
1.9002, 2.0927 

ate.anyは、PSweight()関数によって返されるPSweightオブジェクトである。 ate.anyをprintすると、各処理水準における潜在的な成果の推定平均値のみが得られる。この場合、1.9002 と 2.0927 は、全人口が学歴を持たない場合(Dany = 0)と持たない場合(Dany = 1)の平均対数時間給に相当する。学歴が高いほど平均時給が高くなることがわかる。そのシンプルな画面出力にもかかわらず、ate.anyは、推定傾向スコア(propensity)、推定平均潜在成果(muhat)、推定平均潜在成果の共分散行列(covmu)、bootstrap = TRUEの場合の各 bootstrap標本の推定値(muboot)、アルファベット順のグループラベル(group)、ATTを定義するための指示治療群(trtgrp)の、6要素の一覧から成っている。

PSweight()関数でweightオプションを指定することで、治療集団と重複集団の間の平均的な潜在的成果を同様の方法で推定することができる。 weightが指定されていない場合、PSweight()関数はデフォルトでOWを使用し、最適な内部妥当性を持つ部分母集団をエミュレーションする(Li et al. 2018). ウェイト = "treated "の場合、学歴のある母集団における潜在的成果の平均値を推定して得る。ATTを推定するために、trtgrpを無指定にすると、PSweight()関数はデフォルトで、ATTをdefiingする際に治療変数の最後の値(アルファベット順)を治療グループとみなす(Dany = 1)。 もし、学歴のない母集団での因果関係を推定したい場合は,trtgrp = 0 を指定する必要がある。

ato.any <-  PSweight(ps.formula = ps.any, yname = "wage", data = NCDS)
att.any <-  PSweight(ps.formula = ps.any, yname = "wage", data = NCDS,
                     weight = "treated")
ato.any
Original group value:  0, 1 

Point estimate: 
1.8617, 2.0408 

att.anyの画面出力には、ate.anyやato.anyと比較して、ATT推定値を決定する処置群という要素が含まれている。NCDSデータの分析では、3つの対象集団の平均的な潜在的成果の推定値にわずかな差異が見られるだけであった。 例えば、母集団レベルの効果的な教育介入によって、いずれかの母集団に属するすべての人が学歴を獲得した場合、対数時給の平均値は一貫して高くなるようである。デザインモジュールと同様に、治療効果の加重平均とその分散を推定するためのsummary.PSweight()関数を提供する。

ate.anyおよびato.anyと比較して、att.anyの画面上の出力には、ATT推定値を定義する処置グループであるextra要素が含まれるようになった。NCDSデータの分析では、3つの対象集団における推定平均潜在的転帰の間にわずかな違いしか見られない。平均時間当たりログ賃金は、いずれかの対象集団のすべての個人が、例えば何らかの効果的な集団レベルの教育的介入を通じて学術的資格を獲得した場合、一貫して高くなるように見える。設計モジュールと同様に、概要を示す。PSweight () 関数を使用して、治療効果 (加重平均) とその分散を推定する。既定では、summaryのようになる。PSweight () は、推定された平均の潜在的な結果のtype="DIF")のすべての対ごと対比を示したものであり、したがって、追加の因果的な推定値を対象としている。たとえば、次のコードを使用して、ATEとATO、およびそれらのサンドイッチ標準誤差と95%信頼区間を推定できる。

summary(ate.any, CI  = FALSE)
Closed-form inference: 

Original group value:  0, 1 

Contrast: 
            0 1
Contrast 1 -1 1

           Estimate Std.Error z value  Pr(>|z|)    
Contrast 1 0.192543  0.021122  9.1158 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
summary(ato.any, CI  = FALSE) 
Closed-form inference: 

Original group value:  0, 1 

Contrast: 
            0 1
Contrast 1 -1 1

           Estimate Std.Error z value  Pr(>|z|)    
Contrast 1 0.179129  0.015609  11.476 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

サマリーの戻り値。PSweight () 関数は, ATOの標準誤差がより小さく,関連する信頼区間がよりタイトであり, Liら (2018) の理論的結果と一致することを示す。概要。またPSweight ()関数は、指定された平均の可能性のある結果のコントラストが0であるという弱い帰無仮説のp値を返す。この場合, p値はH_0に対応し,[tex: H{0}=\mu{1}^{n}=\mu_{0}^{n}]はすべて小さく,帰無仮説を棄却する。

傾向スコアモデルを特定することに加えて,各治療群内の潜在的交絡因子の関数として時間給の対数に対するモデルを特定することにより拡張推定量を得る。PSweight ()関数を使用すると、傾向スコアの重み付けと結果のモデリングを組み合わせて、効率性や堅牢性の向上を実現できる。ps.anyで調整されたのと同じ交絡子のセットを使用して、out.formulaを通じて回帰式を指定する。

out.wage <-  wage ~  white + maemp + as.factor(scht) + as.factor(qmab) +
      as.factor(qmab2) + as.factor(qvab) + as.factor(qvab2) + paed_u + maed_u +
      agepa + agema + sib_u + paed_u*agepa + maed_u*agema

処置変数はout.wageには含まれないが、これは、PSweightが各処置群内の個別の潜在的転帰回帰モデルに自動的に適合するためであり、したがって、完全な処置-共変量相互作用が可能となる。連続成果賃金の場合、PSweight () は既定で線形モデルに適合する (family="gaussian") 。結果回帰式を読み込み、augmentation=TRUEを指定すると、セクション2.4で導入された拡張重み付け推定器によって推定された平均の潜在的結果が得られる。

ate.any.aug <-  PSweight(ps.formula = ps.any, yname = "wage", data  = NCDS,
                         weight= "IPW",  augmentation = TRUE, out.formula = out.wage)
ato.any.aug <-  PSweight(ps.formula = ps.any, yname = "wage", data  = NCDS,
                         augmentation = TRUE, out.formula = out.wage)

単純な重み付け推定量と同様に、画面上の出力には、処置群レベルの情報に加えて、それぞれの対象集団における推定された潜在的転帰の平均 (簡略化のために出力を省略) が含まれる。この解析では,検討中の各重みづけ方式に対して,単純重みづけ推定器と拡張重みづけ推定器との間で点推定が実質的に異ならないことを見出した。増強された加重推定値が単純な加重推定値に類似しているという事実は、傾向スコアモデルが大きく誤って特定されていないことの間接的な証拠となる可能性がある (Robins and Rotnitzky 2001; Mercatanti and Li 2014) 。

次に,要約を用いて (加重) 平均治療効果を推定した。PSweight ()関数。この例では、点推定値は、増大重み推定器と単純重み推定器との間で実質的に変化しないが、結果増大は、ATEを推定するための標準誤差を低減するが、ATOを推定するための標準誤差はそれほど大きくない。このような比較結果はMao et al. (2018)のシミュレーション結果と一致する。全体として、考慮された重み付けスキームにかかわらず、平均して学業資格を達成することは、0.05年レベルでないよりも有意に高い時間給につながることがわかった。しかしながら、我々は、研究結果の解釈が0.05より大きいか小さいp値の単一の二分法に依存すべきでないことを認めている。

summary(ate.any.aug, CI=FALSE) 
Closed-form inference: 

Original group value:  0, 1 

Contrast: 
            0 1
Contrast 1 -1 1

           Estimate Std.Error z value  Pr(>|z|)    
Contrast 1 0.186079  0.019842  9.3782 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
summary(ato.any.aug, CI  = FALSE) 
Closed-form inference: 

Original group value:  0, 1 

Contrast: 
            0 1
Contrast 1 -1 1

           Estimate Std.Error z value  Pr(>|z|)    
Contrast 1 0.180004  0.015646  11.505 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

あるいは、標準誤差および信頼区間は、ノンパラメトリックブートストラップを介して推定することができる。たとえば、PSweight () 関数でbootstrap=TRUEを指定し、summaryを使用できる。PSweight () を使用して、平均の潜在的な結果に基づいて、因果的なコントラストのブートストラップベースの推論を行う。デフォルトでは、ブートストラップ複製の数は50に設定されており、PSweight () 関数でR引数を使用して他の値を指定できる。bootstrap=TRUEの場合、PSweight () は50回実行するごとに短いメッセージを出力し、監視を容易する。

ate.any.bs <-  PSweight(ps.formula = ps.any, yname = "wage", data  = NCDS,
                        weight= "IPW",  bootstrap = TRUE)

表示

bootstrap 50 samples

画面上の出力であるate.any.bsは、ate.anyと同じだが、ate.any.bsを要約すると、ブートストラップ標準誤差、 (分位値ベースの) 信頼区間、および関連するp値が返されるようになった。次のコードを使用して、これらの情報を取得する方法を説明する。

summary(ate.any.bs, contrast = rbind(c(-1, 1),c(1, -1)))
Use Bootstrap sample for inference: 

Original group value:  0, 1 

Contrast: 
            0  1
Contrast 1 -1  1
Contrast 2  1 -1

            Estimate Std.Error       lwr      upr  Pr(>|z|)    
Contrast 1  0.192543  0.021948  0.149661  0.22636 < 2.2e-16 ***
Contrast 2 -0.192543  0.021948 -0.226364 -0.14966 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

上記の例では、contrast 引数でデフォルトでない対照を指定できることをさらに説明している。contrast = rbind(c(-1, 1),c(1, -1)) とすることで、Dany = 1 対 Dany = 0 の因果比較とその逆比較を同時に報告することができる。この2つの対比は、同じ因果関係を2つの逆方向から研究しているので、同じ数値が逆符号で返されることが期待される。 ブートストラップ標準誤差はサンドイッチ標準誤差とほぼ同じであるが、ブートストラップ信頼区間は正規近似に依存しないため、点推定値に対して対称でなくなる。

木下優樹菜さん「脳波でADHDが分かった」発言で物議 医師らは全否定

www.msn.com

脳波QEEGをとるのはてんかんをみるためくらいだ。確かに。てんかん型放電はあったのだろうか。

ides.hatenablog.com

ADHDだからといって前頭葉の問題とみるのは早すぎるというか、1/3くらいは小脳に問題を抱えているのにね、とも思う。コンサータが効きました、というなら前頭葉で間違いはないわけだが。
ベタヒスチンの大量投与で小脳タイプのADHDへのアプローチはどうなるんだろうか。アメリカで経口での治験はずいぶん昔に失敗したはずだが、鼻スプレーで治験をするという話が出ている。H3アンタゴニストでどのくらい有効なのかは実際気になるところだ。

精神病性障害患者の犯罪被害リスク要因

www.ncbi.nlm.nih.gov

  • Chapple, B., Chant, D., Nolan, P., Cardy, S., Whiteford, H., & McGrath, J. (2004). Correlates of victimisation amongst people with psychosis. Social Psychiatry and Psychiatric Epidemiology, 39(10), 836–840. https://doi.org/10.1007/s00127-004-0819-4

精神病性障害患者の犯罪被害について調べた研究。962人中172人、17.9 %が犯罪被害を受けている。これは一般人口と比較しても非常に高い。分析結果は2016年に訂正されている(http://link.springer.com/article/10.1007/s00127-015-1147-6)。図表は訂正後のものである。

モデル1は調整なしのオッズ比。モデル2は年齢と性別を調整したオッズ比である。 関連があるのは、年齢(1.62倍)、12か月以内のホームレス体験(3.46倍)、薬物濫用・依存(1.84倍)、12か月以内の逮捕経験(3.68倍)、社会的職業的機能障害(2.73倍)である。社会的職業的機能障害にはSOFASが使用されており、他の項目と合わすためにバイナリで使用されている。高機能のカットオフは70である。SOFASのカットオフについての箇所を引用しよう。

For the purposes of this study, the scores were dichotomised into "slight or no impaired social and occupational functioning” (SOFAS>69) and “poor social and occupational functioning” (SOFAS<70)

SOFAS>69がslight or no impairedで、SOFAS<70がpoorだと書いてあるが、この書き方だと両方70を含んでいるように読める。SOFASは0と1の間でカテゴリカルに分類するので、SOFASは71以上と70以下と分けるのが正しい。この論文筆者は書き間違いをしたか、混乱して理解をしているのかよくわからない(2016年の訂正表では正しく記述されているので現在は正しく理解しているのは間違いない)。

重要なのは、SOFASを2値で使用する際のカットオフの基準は70という点である。。

社会的機能: それは抗うつ薬のエンドポイントになるべきか

www.ncbi.nlm.nih.gov


要約

DSM-IVは社会的機能を測定するための臨床家によって評価される全般的な評価尺度としての社会的・職業的機能尺度(SOFAS)の使用を推奨している。この尺度は、伝統的に抑うつ症状を有する患者を測定する二次的なアウトカム指標として使用される臨床全般印象度(CGI)尺度に類似している。

しかし、 私たちは、抗うつ薬臨床試験におけるエンドポイントとしてこの側面を考えるとき健康に関するQOLが社会的機能の最も適切な指標であると信じている。

健康に関いるQOLは純粋に主観的計測であり、患者による評価がされた評価アンケートは、この文脈の中で最も重要な問題であることが判明している。

この点で、シーハン不安尺度(the Sheehan Disability Scale)は、抗うつ薬臨床試験における社会的機能の最も示唆的な高い全般的な自己報告評価として推奨されている。質問紙のレビューでは、抗うつ薬の試験で社会的機能に関する情報を取得することを目的として3つの最も頻繁に使用される尺度があることを示している。その尺度とは、社会適応尺度セルフレポート版SAS-SR)、社会適応自己評価尺度(SASS)と短文式健康調査(SF-36)である。

しかし、これらの尺度を使用した抗うつ薬プラセボ対照臨床試験の数は、まだ応答性の点で比較するにはあまりにも制限がある。健康に関するQOLは社会的機能は、例えば、身体的健康と精神的健康(認知と感情の両方の問題を含む)など、社会的機能以外の側面が含まれている。

SF-36には、身体的、精神的健康に関連する下位尺度が含まれている。幸福度の観点から測定され、社会的機能は下位尺度になっている。

他のQOLの質問紙である生活の質に関する楽しみと満足感の簡易質問票(Q-LES-Q)のは社会的、精神的、身体的な問題をカバーしている。

最近Q-LES-Qは、60から92の項目からなる総合的な尺度から15項目を含む簡易バージョンへ項目数を減らしている。

追加項目は、全体的な生活の満足度を測定する。

簡便なQ-LES-Qの多くなは社会的機能が含まれており、尺度は、治療の満足に焦点を当てる際には、SF-36やシーハン不安尺度の代替とみなされている。

しかし、これらの質問票を用いた抗うつ薬臨床試験を比較するには数が十分とはいえない。.

このレビューでは、SF-36の下位尺度と抗うつ薬の試験の例では、ほとんどSSRI関連のものだ。.

これらの試験は、抗うつ薬は、6週間の治療期間にわたってプラセボと比較して社会的機能は改善している。エンドポイントスコアはこの時点で国の基準をまだ有意に下回っている。

唯一の12週間の治療は国の基準の限界内に社会的機能のエンドポイント・スコアがある。

再発予防トライアル、もしくはうつ病の再発を防止するための維持療法トライアル、国家の基準を用いた社会的機能尺度の比較では、重要な補足的な指標が治療のために必要がある。

結論として、患者が評価をするQOLの計測の健康に関する一環として社会的機能は、大うつ病の患者における治療の目標を明確にする手助けをするために、抗うつ薬臨床試験における一つのエンドポイントとして必要とされるべきである。

中長期的な試験では、SF-36の下位尺度が、症状志向尺度にサプリメントとして使用されるべきだ。

短期の(6~8週間)期間のトライアルでは、SAS-SRのような他の尺度を使用することは、Q-LES-Qやシーハン不安尺度が考慮されるべきだ。

これらのスケールは、むしろ代替案よりも相互に補完的ととらえるべきである。トライアルの中でこれらのスケールの1つ以上のいずれかを使用する必要があるかもしれないということである。

ベイジアンt検定[R]

www.sumsar.net

bayes.t.test 関数

t.test関数は、4つのバージョンのt検定を実行するために使われる。ここでは、1標本と対の標本の選択肢だけを紹介する。bayes.t.testはt検定のベイジアン的な代替法を実行し、t.test関数と互換性のある関数用法を有している。つまり、t.test(x, conf.int=0.8,mu=1)のようにt-testを実行するだけなら、bayes.を前置するだけですぐに動作するはずである。

bayes.t.testの例として、2002年のNature誌の論文The Value of Bees to the Coffee Harvest (doi:10.1038/417708a, pdf)から引用する。この論文で David W. Roubik は、「熱帯農業の柱である自家受粉のアフリカの低木 Coffea arabicaは昆虫による受粉からは何も得られないと考えられていた」にもかかわらず、ハチはコーヒー収穫にとって重要であると論じている。その論拠となるのが、アフリカミツバチが定着する前の1961年から1980年と、アフリカミツバチがほぼ帰化した1981年から2001年の新世界各国のコーヒー平均収穫量(kg / 10,000m)データである。このデータはミツバチ導入後の収量の増加を示しており、一対の t 検定で分析すると p = 0.04という結果になる。これは、ミツバチがずっと忙しく飛び回っていた旧世界の国々での収穫量の増加と比較すると、一対のt検定ではp =0.232となり、「変化なし」と解釈される。データセットは以下の表とcsvファイル(論文での分析に合わせるため、csvファイルにはカリブ海の島々は含まれていない)にある。

https://www.sumsar.net/files/posts/2014-02-04-bayesian-first-aid-one-sample-t-test/roubik_2002_coffe_yield.csv

このデータセットの分析にt検定を用いることが適切でない理由はいくつかある。t-testでは、各国の地理的な位置は考慮されていないし、どのような「集団」から抽出されたサンプルなのかも明らかではない。また、「相関関係は因果関係を意味しない」という古い決まり文句を唱えたくもなる。確かに、1961年と2001年の間にボリビアで変わったことは、ミツバチの導入以外にたくさんあるはずだ。このような反論を承知で、それでも私は、対応のあるbayes.t.testを紹介するために、この論文を使うつもりである。

まず、論文のオリジナル分析を実行しよう。

d <- read.csv("https://www.sumsar.net/files/posts/2014-02-04-bayesian-first-aid-one-sample-t-test/roubik_2002_coffe_yield.csv")
new_yield_80 <- d$yield_61_to_80[d$world == "new"]
new_yield_01 <- d$yield_81_to_01[d$world == "new"]
t.test(new_yield_01, new_yield_80, paired = TRUE, alternative = "greater")
 Paired t-test

data:  new_yield_01 and new_yield_80
t = 3.1156, df = 12, p-value = 0.004464
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
 628.9914      Inf
sample estimates:
mean of the differences 
               1469.769 

この論文では、片側t検定を使っているので、信頼区間が広くなっている。

では、ベイズ的応急処置の代替案だ。

パッケージのインストール

BayesianFirstAid()JAGSを用いる。rjags()のインストールが必要だ。

install.packages("rjags")

rjags()のインストールがうまくいかない場合にはこちらを参照のこと。
https://gist.github.com/casallas/8411082

Githubからのインストールに必要なdevtools()をインストールしていない場合には、まずインストールが必要である。

install.packages("devtools")

次にGithubからBayesianFirstAid()をインストールする

devtools::install_github("rasmusab/bayesian_first_aid")

BayesianFirstAid()でのベイジアンT検定

library(BayesianFirstAid)
bayes.t.test(new_yield_01, new_yield_80, paired = TRUE, alternative = "less")
 Bayesian estimation supersedes the t test (BEST) - paired samples

data: new_yield_01 and new_yield_80, n = 13

  Estimates [95% credible interval]
mean paired difference: 1402 [376, 2441]
sd of the paired differences: 1688 [895, 2698]

The mean difference is more than 0 by a probability of 0.994 
and less than 0 by a probability of 0.006 

ここでは、平均とSDの推定値が得られ、その横に信頼区間が表示されている。また、平均値の増加率が0より大きい確率は99.4%であることもわかる。古い時代のデータも見てみよう。

old_yield_80 <- d$yield_61_to_80[d$world == "old"]
old_yield_01 <- d$yield_81_to_01[d$world == "old"]
bayes.t.test(old_yield_01, old_yield_80, paired = TRUE)
 Bayesian estimation supersedes the t test (BEST) - paired samples

data: old_yield_01 and old_yield_80, n = 15

  Estimates [95% credible interval]
mean paired difference: 689 [-1214, 2913]
sd of the paired differences: 3521 [1186, 5917]

The mean difference is more than 0 by a probability of 0.77 
and less than 0 by a probability of 0.23 

一般的な関数

すべてのBayesianFirstAid検定には、対応するplot, summary, diagnostics, model.code関数がある以下は、新しい世界のデータを使ったこれらの例である。

plot を使うと、事後分布を直接見ることができる。また、事後分布から抽出されたt分布が散在するヒストグラムの形で、事後予測チェックを得ることができる。モデルの適合度とデータの間に大きな不一致がある場合、先に進む前にもう一度考える必要がある。今のところ、post.pred.チェックは問題なさそうである。   この傾向も収量が増加しているが、推定値の精度は低く、95%信頼区間にはゼロが含まれる。また、新時代の国々と比べて、標準偏差が非常に大きいことも注目すべき点である。

fit <- bayes.t.test(new_yield_01, new_yield_80, paired = TRUE)
plot(fit)

summaryでは、事後分布の詳細な概要が表示される。この概要の小数点以下の桁数は、MCMCで近似しているため、実行ごとに数値が若干飛び交うことに注意。もしこれが気になるようであれば、bayes.t.testを呼び出す際に引数n.iterを使ってMCMCの反復回数を増やすとよいだろう。

summary(fit)
  Data
new_yield_01, n = 13
new_yield_80, n = 13

  Model parameters and generated quantities
mu_diff: the mean pairwise difference between new_yield_01 and new_yield_80 
sigma_diff: the scale of the pairwise difference, a consistent
  estimate of SD when nu is large.
nu: the degrees-of-freedom for the t distribution fitted to the pairwise difference
eff_size: the effect size calculated as (mu_diff - 0) / sigma_diff
diff_pred: predicted distribution for a new datapoint generated
  as the pairwise difference between new_yield_01 and new_yield_80 

  Measures
               mean       sd     HDIlo    HDIup %<comp %>comp
mu_diff    1405.128  525.778   375.531 2469.247  0.006  0.994
sigma_diff 1752.408  472.152   939.557 2739.324  0.000  1.000
nu           29.230   28.502     1.003   85.834  0.000  1.000
eff_size      0.855    0.365     0.143    1.565  0.006  0.994
diff_pred  1410.893 2474.509 -2882.521 5477.636  0.220  0.780

'HDIlo' and 'HDIup' are the limits of a 95% HDI credible interval.
'%<comp' and '%>comp' are the probabilities of the respective parameter being
smaller or larger than 0.

  Quantiles
               q2.5%     q25%   median     q75%   q97.5%
mu_diff      381.766 1068.847 1400.403 1732.760 2478.430
sigma_diff  1004.444 1432.091 1692.692 2000.098 2874.133
nu             2.293    9.439   20.266   39.476  106.273
eff_size       0.180    0.605    0.842    1.081    1.607
diff_pred  -2775.570  189.504 1400.464 2635.695 5612.202

diagnosticsMCMC診断を表示し、プロットする(bayes.binom.testの例と同様)。最後にmodel.codeは、解析を再現する R と JAGS のコードを出力し、R スクリプトに直接コピー&ペーストして、そこから変更することができる。さぁ、始めよう。

model.code(fit)
## Model code for Bayesian estimation supersedes the t test - paired samples ##
require(rjags)

# Setting up the data
x <- new_yield_01 
y <- new_yield_80 
pair_diff <- x - y
comp_mu <-  0 

# The model string written in the JAGS language
model_string <- "model {
  for(i in 1:length(pair_diff)) {
    pair_diff[i] ~ dt( mu_diff , tau_diff , nu )
  }
  diff_pred ~ dt( mu_diff , tau_diff , nu )
  eff_size <- (mu_diff - comp_mu) / sigma_diff
  
  mu_diff ~ dnorm( mean_mu , precision_mu )
  tau_diff <- 1/pow( sigma_diff , 2 )
  sigma_diff ~ dunif( sigma_low , sigma_high )
  # A trick to get an exponentially distributed prior on nu that starts at 1.
  nu <- nuMinusOne + 1 
  nuMinusOne ~ dexp(1/29)
}"

# Setting parameters for the priors that in practice will result
# in flat priors on mu and sigma.
mean_mu = mean(pair_diff, trim=0.2)
precision_mu = 1 / (mad(pair_diff)^2 * 1000000)
sigma_low = mad(pair_diff) / 1000 
sigma_high = mad(pair_diff) * 1000

# Initializing parameters to sensible starting values helps the convergence
# of the MCMC sampling. Here using robust estimates of the mean (trimmed)
# and standard deviation (MAD).
inits_list <- list(
  mu_diff = mean(pair_diff, trim=0.2),
  sigma_diff = mad(pair_diff),
  nuMinusOne = 4)

data_list <- list(
  pair_diff = pair_diff,
  comp_mu = comp_mu,
  mean_mu = mean_mu,
  precision_mu = precision_mu,
  sigma_low = sigma_low,
  sigma_high = sigma_high)

# The parameters to monitor.
params <- c("mu_diff", "sigma_diff", "nu", "eff_size", "diff_pred")

# Running the model
model <- jags.model(textConnection(model_string), data = data_list,
                    inits = inits_list, n.chains = 3, n.adapt=1000)
update(model, 500) # Burning some samples to the MCMC gods....
samples <- coda.samples(model, params, n.iter=10000)

# Inspecting the posterior
plot(samples)
summary(samples) 

例えば、ロバスト性はそれほど大きな問題ではなく、データがt分布ではなく正規分布であると仮定したい場合、このスクリプトを少し修正するだけでよい。モデルコードでは、dtdnormに変更し、nuパラメータを削除し、このようなモデルになる。

model_string <- "model {
  for(i in 1:length(pair_diff)) {
    pair_diff[i] ~ dnorm( mu_diff , tau_diff)
  }
  diff_pred ~ dnorm( mu_diff , tau_diff)
  eff_size <- (mu_diff - comp_mu) / sigma_diff

  mu_diff ~ dnorm( mean_mu , precision_mu )
  tau_diff <- 1/pow( sigma_diff , 2 )
  sigma_diff ~ dunif( sigmaLow , sigmaHigh )
}"

またnuパラメータの監視を外す必要がある。

# The parameters to monitor.
params <- c("mu_diff", "sigma_diff", "eff_size", "diff_pred")

以上。