井出草平の研究ノート

Stataにおける処置効果推定 その1[Stata]

blog.stata.com

処置効果推定量とは、観察データに基づき、ある処置が結果に及ぼす因果的な効果を推定するものである。

本日の投稿では、4つの処置効果推定量について説明する。

  • RA: 回帰調整 Regression adjustment
  • IPW: 逆確率重み付け Inverse probability weighting
  • IPWRA: 逆確率重み付けと回帰調整 Inverse probability weighting with regression adjustment
  • AIPW: 拡張逆確率重み付け Augmented inverse probability weighting

マッチング推定量については、後編に譲る。

処置効果推定量には、魔法のように因果関係を抽出するものはないことに注意しなければならない。観察データの回帰分析同様、因果関係の解釈は、合理的な科学的根拠に基づいて行われなければならない。

はじめに

ここでは、処置と結果について説明する。

処置とは新薬のことであり、その結果は血圧やコレステロールの値である。処置とは、外科的処置であり、その結果は患者の患者の可動性である。処置とは、職業訓練プログラムであり、その成果は雇用や賃金であるかもしれない。処置とは、ある製品の売り上げを伸ばすための広告キャンペーンである可能性もある。

母親の喫煙が出産時の赤ちゃんの体重に影響を与えるかどうかを考えてみよう。このような質問には、観察データを用いてのみ答えることができる。実験では倫理的に問題がある。

観察データの問題点は、被験者が治療を受けるかどうかを選択することである。例えば、ある母親がタバコを吸うか吸わないかを決める。被験者は、処置群と非処置群に自己選択したと言われている。

理想的な世界では、因果関係や治療と結果の関係を検証するために実験を計画します。被験者を治療群と非治療群に無作為に割り当てる。治療法を無作為に割り当てることで、治療と結果の独立性が保証され、分析が非常に簡単になる。

因果推論には、各治療水準でのアウトカムの無条件平均の推定が必要である。我々は、データが観察的か実験的かにかかわらず、治療を受けたという条件付きで各被験者の結果を観察するだけである。実験データでは、治療の無作為割付により、治療が結果に独立であることが保証されているので、観察された治療に対する結果の平均が、関心のある無条件平均を推定する。観察データでは、治療の割り当てプロセスをモデル化する。もし我々のモデルが正しければ,治療割り当てプロセスは,我々のモデル中の共変量に条件づけられたランダムと同じとみなされる.

例を考えてみよう。図1は、Cattaneo (2010)が用いたものと同様の観測データの散布図だ。治療変数は、母親の妊娠中の喫煙状況で、結果は、赤ちゃんの出生体重である。

赤い点は妊娠中にタバコを吸った母親、緑の点は吸わなかった母親を表している。喫煙するかどうかは母親自身が選択したことであり、そのことが分析を複雑にしている。

喫煙した母親としなかった母親の赤ちゃんの平均出生時体重を比較することで、喫煙が出生時体重に及ぼす影響を推定することはできない。なぜだろう?もう一度、グラフを見てみてほしい。喫煙の有無にかかわらず、高齢の母親ほど出生児が重い傾向がある。このデータでは、母親の年齢が高いほど喫煙者である可能性も高い。このように、母親の年齢は治療状況と結果の両方に関係している。では、どうすればいいのだろうか?

RA: 回帰調整推定量 The regression adjustment estimator

回帰調整推定量(RA)は、非ランダムな治療割り付けを考慮するために結果をモデル化する。

私たちは、"喫煙する母親が喫煙しないことを選択していたら、結果はどのように変化しただろうか?" または "喫煙しない母親が喫煙することを選択していたら、結果はどのように変化しただろうか? "と尋ねるかもしれない。もし、これらの反実仮想的な質問に対する答えがわかれば、分析は簡単だ。反実仮想的な結果から観測結果を差し引くだけでいい。

反実仮想的な結果は、治療効果に関する文献では、観察されない潜在的な結果 unobserved potential outcomesと呼ばれている。unobservedという単語は省略されることもある。

我々は、これらの未観測の潜在的な結果の測定値を構築することができ、データは次のようになる。

図2では、観測されたデータを実線で、観測されなかった潜在的な結果を空洞の点で示している。赤い点は、喫煙者がタバコを吸わなかったらどうなるかを示している。緑の点は、非喫煙者が喫煙していた場合の潜在的な結果を示している。

そして、観察されたデータ(実線)を使って、2つの治療群に別々の線形回帰モデルをあてはめることで、観察されない潜在的な結果を推定することができる。

図3では、非喫煙者の回帰線(緑の線)と喫煙者の回帰線(赤の線)を示している。

この2本の線が何を意味しているのかを理解しよう。

図4の左側の緑色の点、Observedと書かれているのは、タバコを吸わない母親の観測値だ。緑の回帰線上のE(y0)と書かれた点は、母親の年齢とタバコを吸わなかったと仮定した場合の赤ちゃんの予想出生体重である。赤の回帰線上のE(y1)という点は、同じ母親が喫煙していた場合の赤ちゃんの出生体重の期待値である。

これらの期待値の差が、治療を受けなかった人たちの共変量特異的治療効果を推定する。

さて、もう一つの反実仮想の問題を見てみよう。

図4の右側の赤い点、Observedと書かれているのは、妊娠中にタバコを吸った母親の観測値です。緑と赤の回帰線上の点は、2つの治療条件下での母親の赤ちゃんの予想出生体重(潜在的な結果)を再び表す。

これらの期待値の差は、治療を受けた人の共変量特異的な治療効果を推定する。

我々は,各被験者について,共変量値に条件付けられた平均治療効果(ATE)を推定することに注意されたい.さらに、この効果は、実際にどの治療を受けたかに関係なく、各被験者について推定する。データ中のすべての被験者に対するこれらの効果の平均が、ATEを推定する。

我々は、図4を使って、受けた治療に関係なく、各被験者が各治療水準で得るであろう結果の予測を動機づけることもできる。この話は、上記の話と類似している。データ中のすべての被験者に対するこれらの予測の平均は、各治療水準に対する潜在的な結果の平均(POMs)を推定する。

推定されたPOMsの差が、上述のATEと同じ推定値であることは、心強いことである。

治療された人のATE(ATET)は、ATEと似ているが、治療群で観測された被験者だけを使用するものである。この治療効果の計算方法は、回帰調整(RA)と呼ばれている。

データセットを開いて、Stataを使って試してみよう。

. webuse cattaneo2.dta, clear
. teffects ra (bweight mage) (mbsmoke), pomeans

我々は、最初の括弧の中に、結果変数とその共変量で結果モデルを指定する。この例では、結果変数は bweight で、唯一の共変量は mage です。

コマンド

  • ra: Regression adjustment
  • pomeans: estimate potential-outcome means

変数

  • bweight: Infant birthweight (grams)
  • mage: Mother's age
  • mbsmoke: 1 if mother smoked

我々は、2番目の括弧の中で、治療モデル(単に治療変数)を指定する。この例では、治療変数 mbsmoke だけを指定する。共変量については、次のセクションで説明する。

このコマンドを入力した結果は以下のとおり。

Iteration 0:   EE criterion =  7.878e-24
Iteration 1:   EE criterion =  8.468e-26

Treatment-effects estimation                    Number of obs      =      4642
Estimator      : regression adjustment
Outcome model  : linear
Treatment model: none
------------------------------------------------------------------------------
             |               Robust
     bweight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
POmeans      |
     mbsmoke |
  nonsmoker  |   3409.435   9.294101   366.84   0.000     3391.219    3427.651
     smoker  |   3132.374   20.61936   151.91   0.000     3091.961    3172.787
------------------------------------------------------------------------------

出力は、すべての母親が喫煙した場合の平均出生体重は3,132g、母親が喫煙しなかった場合の平均出生体重は3,409gとなることを報告する。

POMを差し引くことで、出生体重に対する喫煙のATEを推定することができる。3132.374 - 3409.435 = -277.061. また、teffects raコマンドにateオプションを付けて再実行すれば、標準誤差と信頼区間を得ることができる。

ate: estimate average treatment effect in population; the default

. teffects ra (bweight mage) (mbsmoke), ate

Iteration 0:   EE criterion =  7.878e-24  
Iteration 1:   EE criterion =  5.185e-26  

Treatment-effects estimation                    Number of obs     =      4,642
Estimator      : regression adjustment
Outcome model  : linear
Treatment model: none
----------------------------------------------------------------------------------------
                       |               Robust
               bweight | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-----------------------+----------------------------------------------------------------
ATE                    |
               mbsmoke |
(Smoker vs Nonsmoker)  |  -277.0611   22.62844   -12.24   0.000    -321.4121   -232.7102
-----------------------+----------------------------------------------------------------
POmean                 |
               mbsmoke |
            Nonsmoker  |   3409.435   9.294101   366.84   0.000     3391.219    3427.651
----------------------------------------------------------------------------------------

出力は、我々が手作業で計算したのと同じ、-277.061のATEを報告している。ATEは、各母親が喫煙した場合の出生時体重と、喫煙しない場合の出生時体重の差の平均値である。

ATETはteffects raコマンドでatetオプションをつけて推定することもできるが、ここでは省略する。

atet: estimate average treatment effect on the treated

IPW: 逆確率重み付け Inverse probability weighting

RA 推定法は、非ランダムな治療割り付けを説明するために結果をモデル化する。研究者の中には、処置の割り当てプロセスをモデル化し、結果のモデルを指定しないことを好む人もいます。

私たちのデータでは、喫煙者は非喫煙者よりも高齢である傾向があることが分かっている。また、母親の年齢が出生体重に直接影響するという仮説もある。このことは、図1で観察さ れたが、以下にもう一度示す。

この図から、処理割り当てが母親の年齢に依存していることがわかる。この依存性を調整する方法が欲しいところである。特に、高年齢の緑色と低年齢の赤色のポイントがもっとあればいいと思われる。そうすれば、各群の平均出生時体重は変化するはずだ。それが平均値の差にどのように影響するかはわからないが、差の推定値がより良くなることはわかる。

同様の結果を得るために、低年齢層の喫煙者と高年齢層の非喫煙者のウェイトを高くし、高年齢層の喫煙者と低年齢層の非喫煙者のウェイトを低くするつもりである。

我々は、次のような形式のプロビットまたはロジット・モデルを当てはめる。

Pr(woman smokes) = F(a + b*age)

teffectsは、デフォルトでlogitを使用しますが、ここでは、説明のためにprobitオプションを指定する。

そのモデルを適合させると、データ中の各オブザベーションの予測値 Pr(woman smokes) を得ることができる。これを pi と呼ぶことにする。そして、我々のPOM計算(これは単なる平均計算である)を行う際に、我々は観察結果を重みづけするために、これらの確率を使用する。喫煙者である確率が小さいときに重みが大きくなるように、喫煙者のオブザベーションを1/piで重みづけする。非喫煙者である確率が小さいときに重みが大きくなるように、非喫煙者のオブザベーションを1/(1-pi) で重みづけする。

その結果、図1に代わって次のようなグラフになる。

図5では、大きな円は大きな重みを示している。

このIPW定量を用いてPOMを推定するには、以下のように入力する。

. teffects ipw (bweight) (mbsmoke mage, probit), pomeans

最初の括弧の組は、結果モデルを指定する。これは、この場合、単に結果変数であり、共変量はない。2番目の括弧の組は、処置モデルで、結果変数(mbsmoke)、共変量(この場合、mageだけ)、モデルの種類(probit)が続くものを指定する。

. teffects ipw (bweight) (mbsmoke mage, probit), pomeans

Iteration 0:   EE criterion =  1.007e-21  
Iteration 1:   EE criterion =  1.413e-26  

Treatment-effects estimation                    Number of obs     =      4,642
Estimator      : inverse-probability weights
Outcome model  : weighted mean
Treatment model: probit
------------------------------------------------------------------------------
             |               Robust
     bweight | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
POmeans      |
     mbsmoke |
  Nonsmoker  |   3408.979   9.307838   366.25   0.000     3390.736    3427.222
     Smoker  |   3133.479   20.66762   151.61   0.000     3092.971    3173.986
------------------------------------------------------------------------------

出力結果は、全ての母親が喫煙した場合の平均出生体重は3,133g、母親が誰も喫煙しなかった場合の平均出生体重は3,409gとなることが報告された。

このとき、ATEは-275.5であり、仮に下記のように記した場合、

. teffects ipw (bweight) (mbsmoke mage, probit), ate

Iteration 0:   EE criterion =  2.280e-25  
Iteration 1:   EE criterion =  1.413e-26  

Treatment-effects estimation                    Number of obs     =      4,642
Estimator      : inverse-probability weights
Outcome model  : weighted mean
Treatment model: probit
----------------------------------------------------------------------------------------
                       |               Robust
               bweight | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-----------------------+----------------------------------------------------------------
ATE                    |
               mbsmoke |
(Smoker vs Nonsmoker)  |  -275.5007   22.67711   -12.15   0.000     -319.947   -231.0544
-----------------------+----------------------------------------------------------------
POmean                 |
               mbsmoke |
            Nonsmoker  |   3408.979   9.307838   366.25   0.000     3390.736    3427.222
----------------------------------------------------------------------------------------

標準誤差は22.68、95%信頼区間は [-319.9,231.0] であることを知ることができる。

teffects raと同様に、ATETを求める場合は、teffects ipwコマンドにatetオプションを付けて指定すればいい。

IPWRA: 逆確率重み付けと回帰調整 Inverse probability weighting with regression adjustment

RA推定量は、非ランダムな治療割り付けを考慮した結果をモデル化する。IPW 推定法は、非ランダムな処理割り当てを考慮するために、処理をモデル化するものである。IPWRA推定量では、非ランダムな治療割り付けを考慮するために、結果と処置の両方をモデル化します。

IPWRAは、IPW重みを使用して、その後回帰調整を行うために使用される修正回帰係数を推定する。

結果モデルと処置モデルの共変量は同じである必要はなく、被験者の処置グループの選択に影響する変数が結果に関連する変数と異なることが多いので、そうでないことがよくある。IPWRA推定量にはダブルロバスト特性があり、これは治療モデルまたは結果モデルのどちらか一方(両方ではない)が誤って規定されても、効果の推定値が一致することを意味する。

低出生体重児のデータを使って、より複雑なアウトカムと治療モデルを考えてみよう。

アウトカムモデルには、

  • 1.mage: 母親の年齢
  • 2.prenatal1: 最初の妊娠期間中の出産前訪問を表す指標
  • 3.mmarried: 母親の婚姻状況を表す指標
  • 4.fbaby:第一子であることを表す指標

処置モデルには

  • 1.結果モデルのすべての共変量
  • 2.mage2
  • 3.medu: 母親の教育年数

また、結果モデルと処置モデルの係数を報告するために、aequationsオプションを指定する。

aequations: display auxiliary-equation results

 . teffects ipwra (bweight mage prenatal1 mmarried fbaby)                ///
>                  (mbsmoke mmarried c.mage##c.mage fbaby medu, probit)   ///
>                  , pomeans aequations

Iteration 0:   EE criterion =  9.372e-21  
Iteration 1:   EE criterion =  5.435e-26  

Treatment-effects estimation                    Number of obs     =      4,642
Estimator      : IPW regression adjustment
Outcome model  : linear
Treatment model: probit
-------------------------------------------------------------------------------
              |               Robust
      bweight | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
--------------+----------------------------------------------------------------
POmeans       |
      mbsmoke |
   Nonsmoker  |   3403.336    9.57126   355.58   0.000     3384.576    3422.095
      Smoker  |   3173.369   24.86997   127.60   0.000     3124.624    3222.113
--------------+----------------------------------------------------------------
OME0          |
         mage |   2.893051   2.134788     1.36   0.175    -1.291056    7.077158
    prenatal1 |   67.98549   28.78428     2.36   0.018     11.56933    124.4017
     mmarried |   155.5893   26.46903     5.88   0.000      103.711    207.4677
        fbaby |   -71.9215   20.39317    -3.53   0.000    -111.8914   -31.95162
        _cons |   3194.808   55.04911    58.04   0.000     3086.913    3302.702
--------------+----------------------------------------------------------------
OME1          |
         mage |  -5.068833   5.954425    -0.85   0.395    -16.73929    6.601626
    prenatal1 |   34.76923   43.18534     0.81   0.421    -49.87248    119.4109
     mmarried |   124.0941   40.29775     3.08   0.002     45.11193    203.0762
        fbaby |   39.89692   56.82072     0.70   0.483    -71.46966    151.2635
        _cons |   3175.551   153.8312    20.64   0.000     2874.047    3477.054
--------------+----------------------------------------------------------------
TME1          |
     mmarried |  -.6484821   .0554173   -11.70   0.000     -.757098   -.5398663
         mage |   .1744327   .0363718     4.80   0.000     .1031452    .2457202
              |
c.mage#c.mage |  -.0032559   .0006678    -4.88   0.000    -.0045647   -.0019471
              |
        fbaby |  -.2175962   .0495604    -4.39   0.000    -.3147328   -.1204595
         medu |  -.0863631   .0100148    -8.62   0.000    -.1059917   -.0667345
        _cons |  -1.558255   .4639691    -3.36   0.001    -2.467618   -.6488926
-------------------------------------------------------------------------------

出力のPOmeansセクションは、2つの処理グループのPOMを表示します。ATEは、3173.369 - 3403.336 = -229.967と計算され。

OME0セクションとOME1セクションは、それぞれ未処理群と処理群のRA係数を表示する。

出力のTME1 セクションは、probit 処置モデルの係数を表示する。

前の2つのケースと同様に、もし、標準誤差などを含むATEが欲しければ、ateオプションを指定する。ATETが欲しい場合は、atetオプションを指定する。

AIPW: 拡張逆確率重み付け Augmented inverse probability weighting

IPWRA推定量では、非ランダムな治療割り付けを考慮し、結果と治療の両方をモデル化する。AIPW推定量も同様である。

AIPW推定量は、IPW定量にバイアス補正項を追加する。治療モデルが正しく規定されていれば、バイアス補正項は0であり、モデルはIPW定量に縮小される。もし治療モデルが誤って規定されていても、結果モデルが正しく規定されていれば、バイアス補正項が推定量を修正する。このように、バイアス補正項はAIPW推定量にIPWRA推定量と同じダブルロバスト特性を与える。

AIPW推定量の構文と出力は、IPWRA推定量とほぼ同じである。

 teffects aipw (bweight mage prenatal1 mmarried fbaby)                 ///
>               (mbsmoke mmarried c.mage##c.mage fbaby medu, probit)    ///
>               , pomeans aequations

Iteration 0:   EE criterion =  4.632e-21  
Iteration 1:   EE criterion =  5.391e-26  

Treatment-effects estimation                    Number of obs     =      4,642
Estimator      : augmented IPW
Outcome model  : linear by ML
Treatment model: probit
-------------------------------------------------------------------------------
              |               Robust
      bweight | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
--------------+----------------------------------------------------------------
POmeans       |
      mbsmoke |
   Nonsmoker  |   3403.355   9.568472   355.68   0.000     3384.601    3422.109
      Smoker  |   3172.366   24.42456   129.88   0.000     3124.495    3220.237
--------------+----------------------------------------------------------------
OME0          |
         mage |   2.546828   2.084324     1.22   0.222    -1.538373    6.632028
    prenatal1 |   64.40859   27.52699     2.34   0.019     10.45669    118.3605
     mmarried |   160.9513    26.6162     6.05   0.000     108.7845    213.1181
        fbaby |   -71.3286   19.64701    -3.63   0.000     -109.836   -32.82117
        _cons |   3202.746   54.01082    59.30   0.000     3096.886    3308.605
--------------+----------------------------------------------------------------
OME1          |
         mage |  -7.370881    4.21817    -1.75   0.081    -15.63834    .8965804
    prenatal1 |   25.11133   40.37541     0.62   0.534    -54.02302    104.2457
     mmarried |   133.6617   40.86443     3.27   0.001      53.5689    213.7545
        fbaby |   41.43991   39.70712     1.04   0.297    -36.38461    119.2644
        _cons |   3227.169   104.4059    30.91   0.000     3022.537    3431.801
--------------+----------------------------------------------------------------
TME1          |
     mmarried |  -.6484821   .0554173   -11.70   0.000     -.757098   -.5398663
         mage |   .1744327   .0363718     4.80   0.000     .1031452    .2457202
              |
c.mage#c.mage |  -.0032559   .0006678    -4.88   0.000    -.0045647   -.0019471
              |
        fbaby |  -.2175962   .0495604    -4.39   0.000    -.3147328   -.1204595
         medu |  -.0863631   .0100148    -8.62   0.000    -.1059917   -.0667345
        _cons |  -1.558255   .4639691    -3.36   0.001    -2.467618   -.6488926
-------------------------------------------------------------------------------

ATEは3172.366 - 3403.355 = -230.989となる。

最終的な考察

上記の例では、出生体重という連続的なアウトカムを使用した。teffectsは、二値、計数、非負の連続的なアウトカムでも使用することができる。

また、推定量では、複数の治療カテゴリーが可能である。

Stata 13のtreatment-effects機能については、マニュアル全体が費やされており、基本的な紹介、高度な議論、そして実例が含まれている。さらに詳しく知りたい方は、Stataのウェブサイトから[TE] Treatment-effects Reference Manualをダウンロードすることができる。 http://www.stata.com/manuals13/te/

参考文献

Cattaneo, M. D. 2010. Cattaneo, M. D. 2010: Efficient semiparametric estimation of multi-valued treatment effects under ignorability. Journal of Econometrics 155: 138-154.

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

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つ以上のいずれかを使用する必要があるかもしれないということである。