日誌

Stataメモ

Stataメモ目次

参照に便利なように,目次を作りました

researchmapのバージョンアップによるリンク切れを修正 2020223

【データハンドリング】


【解析】


【図・グラフ】

0

複数の回帰分析の結果をtableで報告する形式に出力する

回帰分析で説明変数が1つだけの2変量モデルから,調整する変数
を段階的に加えていき,複数のモデルを並べてtableで報告することが
よく行われます。普通にやったら手動で編集して作表しなければ
なりませんが,ユーザー作成コマンドのestoutを使うととても
楽にできます。

まずパッケージestoutをインストールします。まだ入ってない場合,
ネットにつながっている状態で次のように打ちます。

ssc install estout


autoデータセットを使います。
rep78*foreignのクロス表では,0名のセルが発生するので,
単純にするためにrep78の3未満をdropしておきます。

sysuse auto
drop if rep78<3


(1)回帰係数のみを出力

まずはアウトカムをforeign,説明変数をmpgだけにしたロジスティック回帰

logistic foreign mpg

--------------------------------------------------------------------
foreign | Odds Ratio   Std. Err.    z    P>|z|  [95% Conf. Interval]
--------+-----------------------------------------------------------
    mpg |   1.160921   .0614601   2.82   0.005    1.0465    1.287852
  _cons |    .019403   .0235867  -3.24   0.001  .0017911    .2101887
--------------------------------------------------------------------

ここで,次のコマンドを実行すると,回帰係数だけが出てきます。
ただし,オッズ比ではありません。

estout

-------------------------
                        b
-------------------------
foreign                 
mpg              .1492138
_cons           -3.942329
-------------------------

念のため,オッズ比変換前の回帰係数も出してみて確認

logit foreign mpg

--------------------------------------------------------------------
foreign |      Coef.   Std. Err.    z    P>|z|  [95% Conf. Interval]
--------+-----------------------------------------------------------
    mpg |   .1492138   .0529408   2.82   0.005  .0454517    .2529758
  _cons |  -3.942329   1.215624  -3.24   0.001 -6.324908    -1.55975
--------------------------------------------------------------------


(2)複数の回帰分析の結果を出力

さて,それでは,1つずつ説明変数を加えて行った,計3つのモデルを
作成します。出力結果を出さないように,quietlyを追加しています。

eststo:quietly logistic foreign mpg
eststo:quietly logistic foreign mpg length
eststo:quietly logistic foreign mpg length i.rep78, allbaselevels

これで3つのモデルが格納されたので確認してみます。

eststo dir

        name | command      depvar       npar  title
-------------+-----------------------------------------
        est1 | logistic     foreign         2 
        est2 | logistic     foreign         3 
        est3 | logistic     foreign         6 
-------------------------------------------------------


3つのモデルすべてを出力します。*はすべてのモデルを指定する記号です。

estout *

                     est1         est2         est3
                        b            b            b
---------------------------------------------------
foreign                                           
mpg              .1492138    -.0900729    -.1824287
length                       -.0963129    -.1063041
3.rep78                                           0
4.rep78                                    2.405568
5.rep78                                    3.729435
_cons           -3.942329     18.82339     21.00963

一列目が説明変数で,モデルごとに列が変わり,説明変数が増えていく
様が分かります。


(3)出力結果を整える

出したい形式は,コマンドが複雑なので,ちょっとずつ
オプションを追加して説明していきます。

・回帰係数と95%信頼区間を表示し,切片を削除

estout *, cell(b ci)  drop(_cons)

---------------------------------------------------
                   est1                est2         est3
                 b/ci95              b/ci95       b/ci95
---------------------------------------------------
foreign                                          
mpg            .1492138          -.0900729        -.1824287
          .0454517,.2529758 -.2432507,.0631049 -.3959393,.0310819
length                           -.0963129         -.1063041
                            -.1529248,-.0397009 -.1725313,-.0400769
3.rep78                                                  0
                                                       0,0
4.rep78                                           2.405568
                                              .3872429,4.423893
5.rep78                                           3.729435
                                              1.175191,6.283679
---------------------------------------------------

それぞれの回帰係数の下の行に信頼区間が追加されています

・オッズ比とオッズ比の95%信頼区間

estout *, cells(b ci) drop(_cons) eform

                     est1            est2             est3
                   b/ci95          b/ci95           b/ci95
----------------------------------------------------------
foreign                                                 
mpg              1.160921        .9138646         .8332441
             1.0465,1.287852 .7840749,1.065139 .6730476,1.03157
length                           .9081798         .8991512
                             .8581942,.9610768 .8415319,.9607155
3.rep78                                                 1
                                                      1,1
4.rep78                                          11.08472
                                             1.472914,83.42037
5.rep78                                          41.65555
                                             3.23876,535.756


・小数点第1位まで表示し,切片を削除

estout *, cells(b(fmt(1)) ci(fmt(1))) drop(_cons) eform

---------------------------------------------------
                     est1         est2         est3
                   b/ci95       b/ci95       b/ci95
---------------------------------------------------
foreign                                           
mpg                   1.2          0.9          0.8
                  1.0,1.3      0.8,1.1      0.7,1.0
length                             0.9          0.9
                               0.9,1.0      0.8,1.0
3.rep78                                         1.0
                                            1.0,1.0
4.rep78                                        11.1
                                           1.5,83.4
5.rep78                                        41.7
                                          3.2,535.8


・95%信頼区間を横に配置
bからciの部分を"  "で囲みます

estout *, cells("b(fmt(1)) ci(fmt(1))") drop(_cons) eform

         est1                est2                est3         
            b       ci95        b       ci95        b       ci95
----------------------------------------------------------------
foreign                                                       
mpg       1.2    1.0,1.3      0.9    0.8,1.1      0.8    0.7,1.0
length                        0.9    0.9,1.0      0.9    0.8,1.0
3.rep78                                           1.0    1.0,1.0
4.rep78                                          11.1   1.5,83.4
5.rep78                                          41.7  3.2,535.8


・95%信頼区間を[ ]で閉じる
おそらくこれが一番有用な形式と思います。ciの所にparを追加します

estout *, cells("b(fmt(1)) ci(par fmt(1))") eform drop(_cons)

         est1               est2              est3           
            b       ci95       b       ci95      b         ci95
---------------------------------------------------------------
foreign                                                      
mpg       1.2  [1.0,1.3]     0.9  [0.8,1.1]    0.8    [0.7,1.0]
length                       0.9  [0.9,1.0]    0.9    [0.8,1.0]
3.rep78                                        1.0    [1.0,1.0]
4.rep78                                       11.1   [1.5,83.4]
5.rep78                                       41.7  [3.2,535.8]



・95%信頼区間を( )で閉じる
parの所に( , )を追加します。

estout *, cells("b(fmt(1)) ci(par(( , )) fmt(1))") eform drop(_cons)

         est1               est2              est3           
            b       ci95       b       ci95      b         ci95
---------------------------------------------------------------
foreign                                                      
mpg       1.2  (1.0,1.3)     0.9  (0.8,1.1)    0.8    (0.7,1.0)
length                       0.9  (0.9,1.0)    0.9    (0.8,1.0)
3.rep78                                        1.0    (1.0,1.0)
4.rep78                                       11.1   (1.5,83.4)
5.rep78                                       41.7  (3.2,535.8)
0

介入研究の解析での線形混合モデル(1)

本記事の解説では,Stata 13.1を使います。

ここで解説するRCTなどの介入研究で使われる解析法は,
線形混合モデル(linear mixed model),マルチレベル(multilevel),
階層線型モデル(hierarchical linear model)など様々な呼び方が
ありますが,ここでは介入研究の文脈でよく用いられる線形混合モデル
の呼び方にします(「線型」という表記もあります)。

従来は,アウトカムが連続量であるRCTなどの3時点以上の全体の解析は
反復測定分散分析(repeated ANOVA)で行われていました。しかし近年では,
これが線形混合モデルになることが多いので,こうしたタイプの
解析を行う研究者にとって必須の知識になりつつあると思います。

特に,追跡時に必ずと言っていいほど発生する欠損について,
データがそろっている対象者だけに限定するとか,前の時点の値を代入する
(LOCF)とかをしなくても,問題なく解析できる点が強みです。

LOCFについては↓参照
追跡データの欠損について前の時点の値を代入する(LOCF)

いろいろ勉強しないとわからないことが多いのですが,まず解析
して体験することでわかることも多いと思いますので,
これまで解説に使ってきたBtheBデータを使って試してみます。

ちなみにBtheBデータセットは,実際の研究のデータで,↓のように
論文も出てるので,解釈の参考にもなります。


BtheBデータの使い方は下記のページ↓参照

EZRでパッケージ中のRデータセットをStataデータセットに変換する

BtheBデータにある欠損値の状況は下記のページ↓で解説してます
サンプルデータで現実的な欠損値が含まれているものは少ない気がする
ので,これがこのデータセットの強みです。

データの欠損を数える

線型混合モデルで縦断データを解析する際の準備は,下記のページの作業↓
が必要になります

データセットをwideからlongへ変換する

さて,これで準備が整いましたので,いよいよ解析に進みます。
基本的な解析自体は非常に簡単で,抑うつ症状(bdi)について
群(treatment)×時間(month)の交互作用の検証という主目的は以下の
1行だけです。以下はランダム切片だけですが,Stataだとこんなシンプルに
書けるのが素敵です。

mixed bdi treatment##month || id:

##の記号で説明変数をつなぐだけで,
それぞれの主効果も自動で含めてくれます。Stataではバージョン12まではxtmixedという
コマンドでしたが,13からmixedになっており,徐々に進化している感があります。

詳しくはhelp mixedと打つといろいろわかります。

出力で最も関心のある交互作用を確認するには↓のtreatment#monthの所を見ます。

※2016/01/31 修正 オムニバスの交互作用は最近ではあまり使われないです。詳しくは
介入研究の解析での線形混合モデル(3)介入効果はどこを見るか 参照

----------------------------------------------------------------------------------------------------
以下は記事修正前に書いていた,オムニバス交互作用の検定をしたことを前提に示した
結果となります。
 
----------------------------------------------------------------------------------------------------------------
               bdi |      Coef.        Std. Err.      z      P>|z|     [95% Conf. Interval]
----------------+----------------------------------------------------------------------------------------------
      treatment |
          BtheB  |  -4.755128   2.211978    -2.15   0.032    -9.090525   -.4197319
                     |
          month |
                 3  |  -1.606839   1.160547    -1.38   0.166     -3.88147     .667791
                 5  |  -3.214118   1.259442    -2.55   0.011     -5.68258   -.7456561
                 8  |   -5.94658    1.328405    -4.48   0.000    -8.550205   -3.342955
                     |
treatment#month |
       BtheB#3  |   .5951418    1.62632      0.37   0.714    -2.592386     3.78267
       BtheB#5  |   1.335008   1.774927     0.75   0.452    -2.143785      4.8138
       BtheB#8  |   3.325952   1.847011     1.80   0.072    -.2941233    6.946027
                      |
            _cons |   19.46667   1.619558    12.02   0.000     16.29239    22.64094
--------------------------------------------------------------------------------------------------------------------

デフォルトだとmonthがカテゴリ変数として認識されていました。
これだけでは,論文に記載することの多い交互作用項全体のp値はわかりません。
これを出すためには,以下のコマンドを実行します。

contrast treatment#month

---------------------------------------------------------------------
                |         df        chi2     P>chi2
----------------+-------------------------------------------
bdi             |
treatment#month |          3        3.46     0.3266
---------------------------------------------------------------------


これによって,群×時間の交互作用はp=0.3266で有意ではないという
結果が分かりました。

このカイ二乗値をF値にしたいときは,↓を参考に

How can I test simple effects in repeated measures models? (Stata 12)

カイ二乗値を自由度で割る計算をするだけです。すなわち

display 3.46/3

というように上のカイ二乗値を自由度で割るコマンドを書くと,

1.1533333

という結果が出てきます。これがF値になります。

ちなみに,Stata公式がカイ二乗値とF値の関係について解説しています



入門知識はこちらにまとめていますのでご参考まで

第2回心理・医学系研究者のためのデータ解析環境Rによる統計学の研究会(土屋)-公開用

関連記事:
介入研究の解析での線形混合モデル(2):推定結果のグラフ化


このトピックスは適宜修正・追加予定です。

【履歴】
2014/06/02 記事作成
2014/07/21 一部削除→交互作用の部分
2016/01/31 修正。出力のポイントについて説明と参照リンク追加
2016/02/27 前回修正した部分の表現を修正
0

割合と信頼区間のグラフを描く

割合の95%信頼区間を算出する

にて,割合と信頼区間の算出を紹介しました。
これらをStataでをグラフにするには,
割合と信頼区間それぞれが変数となった新たなデータセットが必要です。
以下,バージョンはStata14.1を用いました。

上記の解析で得られた結果出力をexcel等で加工して取り出して
作ることもできますが,手間がかかります。
そこで,何とかStata上だけで処理する方法を考えてみました。


1. 割合と信頼区間の情報のみを抽出

例えば,rep78=3だけ取り出してみてみます。

ci proportions foreign if rep78==3
 
                                         -- Binomial Exact --
Variable | Obs  Proportion    Std. Err.  [95% Conf. Interval]
-------------+---------------------------------------------------
 foreign |  30          .1    .0547723   .0211171    .2652885


この解析が行われると,いくつかの情報がメモリに格納されるようです。

return list

と打つと出てきます。

scalars:
                  r(N) =  30
         r(proportion) =  .1
                 r(se) =  .0547722557505166
                 r(lb) =  .0211171370297225
                 r(ub) =  .2652884504742086
              r(level) =  95

macros:
             r(citype) : "exact"


この内,割合と信頼区間の情報である

         r(proportion) =  .1
                 r(lb) =  .0211171370297225
                 r(ub) =  .2652884504742086

さえ選んで出力できれば便利です。以下のようにしてできます。
出力したい値をコンマ(,)で区切って並べると,スペースで
区切って出力されます

display r(proportion),r(lb),r(ub)

.1 .02111714 .26528845

あとは,これをrep78の値ごとに繰り返せば
便利です。そんな時はloopを使います。くわしくは
たくさんの変数に同じ処理を繰り返したいとき

forvalues i=3/5 {
 quietly ci proportions foreign if rep78==`i'
 display `i',r(proportion),r(lb),r(ub)
 }


rep78が3-5まであるので,それぞれごとに
foreignの割合と信頼区間を算出し,算出結果の
割合,95%信頼区間の下限,上限を,rep78の値ごと(i)
にdisplayで選んで出力。なお,ciの出力はquietlyで抑制

出力は以下の通り

3 .1 .02111714 .26528845
4 .5 .26019058 .73980942
5 .81818182 .48224415 .9771688

あとはこれをStataデータセットに読み込めれば,
グラフ化に便利です。


2. 新たなデータセット作成

ここでいったんデータセットをclearします。

clear

1.で取り出した情報を結果画面からコピーして
以下のようにエディタに書きます。

input rep78 proportion lower upper
/*↓に貼り付け*/
3 .1 .02111714 .26528845
4 .5 .26019058 .73980942
5 .81818182 .48224415 .9771688
end


inputの行は変数の名前です。そして最後にendを置きます。

次に,100%換算にするために,割合と信頼区間
それぞれに100を掛けます。

replace proportion=proportion*100
replace lower=lower*100
replace upper=upper*100


そして,新たなデータセットを確認します。

list _all

     | rep78   propor~n      lower      upper |
     |----------------------------------------|
  1. |     3         10   2.111714   26.52884 |
  2. |     4         50   26.01906   73.98094 |
  3. |     5   81.81818   48.22441   97.71688 |

これで準備完了です。


3. 割合と信頼区間のグラフ作成(縦書き)

ここまでできたらあとは簡単です。
rspikeで信頼区間,scatterで割合をそれぞれ描画します。

twoway rspike upper lower rep78 || scatter proportion rep78




見栄えが悪いので,オプションを色々追加してグラフを表示します。

twoway rspike upper lower rep78 || scatter proportion rep78,  ///
    msymbol(square) msize(medium) xlabel(3(1)5)  ///
    legend(pos(5) ring(0) col(1) label(1 "95%CI"))  plotregion(margin(large))


msymbolで点推定値を四角形に指定し,msizeで大きさを変更,
xlabelで,rep78が3~5まで1ずつ表示されるようにし,
legendで凡例の位置を調整,ラベルも変更,plotregionで左右の空きを増やしました。




rspikeをrcapに変えると,信頼区間の端に棒がつきます。

twoway rcap upper lower rep78 || scatter proportion rep78,  ///
   msymbol(square) msize(medium)  xlabel(3(1)5)  ///
   legend(pos(5) ring(0) col(1) label(1 "95%CI")) plotregion(margin(large))





4. 割合と信頼区間のグラフ作成(横書き)

横向きに描きたいときは,rcapの方にhorizontalオプションを付け,
scatterの変数の順序を逆にします。

twoway rspike upper lower rep78, horizontal || scatter rep78 proportion

0

割合の95%信頼区間を算出する

stata14.1を使います。使用するコマンドci proportionsは
バージョン14.1から追加されたみたいです。
以前はciでbinomialオプションで扱っていたようです。

autoデータを使います。

複数変数の層別平均値を変数として追加し各層の最初の1行だけ表示

の記事でみたように,rep78とforeignのクロス表では0セルが
発生するので,rep78の1と2は削除しておきます。

sysuse auto
drop if rep78==1 | rep78==2



1.割合の95%信頼区間

まず変数foreignを使います。
0:Domestic, 1:Foreignという値なので,
1:Foreignの割合を算出するということになります。

まずは確認で度数分布

 tab foreign


   Car type |      Freq.     Percent        Cum.
------------+-----------------------------------
   Domestic |         42       65.63       65.63
    Foreign |         22       34.38      100.00
------------+-----------------------------------
      Total |         64      100.00

rep78の1と2をdropしているので人数は減っています。
次に割合と信頼区間

ci proportions foreign

                                          -- Binomial Exact --
 Variable |  Obs  Proportion   Std. Err. [95% Conf. Interval]
----------+-------------------------------------------------
  foreign |   64      .34375   .0593699  .2294632    .473023


proportionが0.34なので,34%です。信頼区間は0.229, 0.473
defaultの信頼区間はexact(Clopper-Pearson)で,
オプションでwald, wilson, agresti, jeffreysに対応しています


2.層別に割合の95%信頼区間

rep78の値ごとに出したい場合は

bysort rep78: ci proportions foreign

-> rep78 = 3

                                           -- Binomial Exact --
   Variable |  Obs  Proportion    Std. Err.  [95% Conf. Interval]
------------+---------------------------------------------------
    foreign |   30          .1    .0547723   .0211171    .2652885

----------------------------------------------------------------
-> rep78 = 4

                                            -- Binomial Exact --
  Variable |  Obs  Proportion    Std. Err.  [95% Conf. Interval]
-----------+----------------------------------------------------
   foreign |   18          .5    .1178511   .2601906    .7398094

-----------------------------------------------------------------
-> rep78 = 5

                                            -- Binomial Exact --
  Variable |  Obs  Proportion    Std. Err.  [95% Conf. Interval]
-----------+----------------------------------------------------
   foreign |   11    .8181818    .1162913   .4822441    .9771688

-----------------------------------------------------------------
-> rep78 = .

                                           -- Binomial Exact -    Variable |  Obs  Proportion    Std. Err.  [95% Conf. Interval]
---------+----------------------------------------------------
 foreign |    5          .2    .1788854   .0050508    .7164179


最後にrep78が欠損の場合も出力されています
0

介入研究の解析での線形混合モデル(3)介入効果はどこを見るか

線形混合モデルを使用した際の、治療効果や介入効果の求め方について、
以下の本,

経時的繰り返し測定デザイン : ─治療効果を評価する混合効果モデルとその周辺─ (医学統計学シリーズ)
丹後 俊郎
朝倉書店(2015/09/19)
値段:¥ 4,860



が非常に詳細に解説していて,とても理解が進んでよかったので、
本の中ではSASで示されていたコード部分の例をStataで再現してみます。
Stataのバージョンは13.1です。

BtheBデータのStataへの読み込み方はたとえば↓の過去記事参照


読み込んだらまず,後にデータセット構造をwideからlongに変換する必要が
あるので,Stataのルール上,変数名を次のように変換します。

rename (bdi_pre bdi_2m bdi_3m bdi_5m bdi_8m) (bdi0 bdi2 bdi3 bdi5 bdi8)

次にIDがないので作成します。

egen id=seq()

最初の5名のデータを確認します。

list in 1/5

     | drug   length   treatm~t   bdi0   bdi2   bdi3   bdi5   bdi8   id |
     |------------------------------------------------------------------|
  1. |   No      >6m        TAU     29      2      2      .      .    1 |
  2. |  Yes      >6m      BtheB     32     16     24     17     20    2 |
  3. |  Yes      <6m        TAU     25     20      .      .      .    3 |
  4. |   No      >6m      BtheB     21     17     16     10      9    4 |
  5. |  Yes      >6m      BtheB     26     23      .      .      .    5 |


(1)アウトカムに1つでも欠損のある者を除いた完全ケースデータを作成(p74のModel II)

keep if bdi0~=. & bdi2~=. & bdi3~=. & bdi5~=. & bdi8~=.
list in 1/5

 
    | drug   length   treatm~t   bdi0   bdi2   bdi3   bdi5   bdi8   id |
     |------------------------------------------------------------------|
  1. |  Yes      >6m      BtheB     32     16     24     17     20    2 |
  2. |   No      >6m      BtheB     21     17     16     10      9    4 |
  3. |  Yes      <6m      BtheB      7      0      0      0      0    6 |
  4. |  Yes      <6m        TAU     17      7      7      3      7    7 |
  5. |   No      >6m        TAU     20     20     21     19     13    8 |


idの1,3,5番は欠損値があったので,データセットから除かれています。
次に群ごとの人数を確認します。

tab treatment

  treatment |      Freq.     Percent        Cum.
------------+-----------------------------------
        TAU |         25       48.08       48.08
      BtheB |         27       51.92      100.00
------------+-----------------------------------
      Total |         52      100.00

p68にある群ごとの人数と一致します。


(2)データセットをwideからlongに変換

reshape long bdi, i(id) j(month)
list in 1/10

     | id   month   drug   length   treatm~t   bdi |
     |---------------------------------------------|
  1. |  2       0    Yes      >6m      BtheB    32 |
  2. |  2       2    Yes      >6m      BtheB    16 |
  3. |  2       3    Yes      >6m      BtheB    24 |
  4. |  2       5    Yes      >6m      BtheB    17 |
  5. |  2       8    Yes      >6m      BtheB    20 |
     |---------------------------------------------|
  6. |  4       0     No      >6m      BtheB    21 |
  7. |  4       2     No      >6m      BtheB    17 |
  8. |  4       3     No      >6m      BtheB    16 |
  9. |  4       5     No      >6m      BtheB    10 |
 10. |  4       8     No      >6m      BtheB     9 |


(3)線形混合モデル(p74のModel II) ※2016/01/31 本の解析と合うように修正

本ではベースライン以外の測定時を1とした変数postを
作成しているので同様に作ります。

gen post=0
replace post=1 if month~=0

本と同様に,treatmentを連続変数として扱うため以下の
様に値を0,1に変換します。

replace treatment=treatment-1 //TAU=1,BthB=2なので,0-1に変換

次が線形混合モデル本体です。 

mixed bdi drug length c.treatment##month || id:post , reml  cov(unstruct) allbase

allbaseオプションは,referenceカテゴリーが何か明示するためにつけてます。
treatmentは連続変数として扱うため,c.をつけています。
カテゴリ変数のreferenceは,Stataではデフォルトが最初のカテゴリなので,SAS
のように変更する必要はありません。

ランダム効果部分の出力

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Unstructured             |
                   var(post) |    61.8672   18.85135      34.04816    112.4158
                  var(_cons) |   54.45944    16.9922      29.54506    100.3833
             cov(post,_cons) |  -28.03576   15.18228     -57.79248    1.720963
-----------------------------+------------------------------------------------
               var(Residual) |   24.56581   2.836617      19.59038    30.80488
------------------------------------------------------------------------------
LR test vs. linear model: chi2(3) = 131.47                Prob > chi2 = 0.0000

母数効果の部分の出力

-----------------------------------------------------------------------------------
              bdi |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
------------------+----------------------------------------------------------------
             drug |   1.123422   2.027048     0.55   0.579    -2.849519    5.096363
           length |   7.225293   1.988499     3.63   0.000     3.327907    11.12268
        treatment |  -1.816101   2.534019    -0.72   0.474    -6.782687    3.150485
                  |
            month |
               0  |          0  (base)
               2  |   3.068148   4.676631     0.66   0.512    -6.097881    12.23418
               3  |  -.8933333   4.676631    -0.19   0.849    -10.05936    8.272696
               5  |  -3.881481   4.676631    -0.83   0.407    -13.04751    5.284547
               8  |  -7.891852   4.676631    -1.69   0.092    -17.05788    1.274177
                  |
month#c.treatment |
               0  |          0  (base)
               2  |  -7.108148   2.924213    -2.43   0.015     -12.8395   -1.376796
               3  |  -5.386667   2.924213    -1.84   0.065    -11.11802    .3446852
               5  |  -4.318519   2.924213    -1.48   0.140    -10.04987    1.412833
               8  |  -2.628148   2.924213    -0.90   0.369      -8.3595    3.103204
                  |
            _cons |    12.6037   5.643376     2.23   0.026      1.54289    23.66452
-----------------------------------------------------------------------------------



この各交互作用項の推定値が,TAU群に比べたBtheB群のベースラインからの変化について
の各測定時点での群間差になり,つまり効果量(非標準化)
ということになります。
つまり,大抵の論文でのメインで報告する値です。どの測定時を主要評価時点
とするかは事前に設定しておく必要があります。

TAU群に比べたBtheB群のベースラインからの変化(介入期間を通じた平均):meanCFB

contrast {c.treatment#month -1 0.25 0.25 0.25 0.25}

-------------------------------------------------------------------
                  |   Contrast   Std. Err.     [95% Conf. Interval]
------------------+------------------------------------------------
bdi               |
month#c.treatment |
             (1)  |   -4.86037   2.670517     -10.09449    .3737474
-------------------------------------------------------------------




--------------------------------------------------------------------------------------------------------
以下はtreatmentがカテゴリ変数だった場合の結果(2016/01/31以前の記事で説明
していた内容を要約)

mixed bdi drug length treatment##month || id:post , reml cov(unstruct) allbase

ランダム効果部分の出力

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Unstructured             |
                   var(post) |    61.8672   18.85135      34.04816    112.4158
                  var(_cons) |   54.45944    16.9922      29.54506    100.3833
             cov(post,_cons) |  -28.03576   15.18228     -57.79248    1.720962
-----------------------------+------------------------------------------------
               var(Residual) |   24.56581   2.836617      19.59038    30.80488
------------------------------------------------------------------------------


母数効果の部分の出力

            bdi |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
----------------+----------------------------------------------------------------
           drug |   1.123422   2.027048     0.55   0.579    -2.849519    5.096363
         length |   7.225293   1.988499     3.63   0.000     3.327907    11.12268
                |
      treatment |
           TAU  |          0  (base)
         BtheB  |  -1.816101   2.534019    -0.72   0.474    -6.782687    3.150485
                |
          month |
             0  |          0  (base)
             2  |      -4.04    2.10712    -1.92   0.055    -8.169879    .0898786
             3  |      -6.28    2.10712    -2.98   0.003    -10.40988   -2.150121
             5  |       -8.2    2.10712    -3.89   0.000    -12.32988   -4.070121
             8  |     -10.52    2.10712    -4.99   0.000    -14.64988   -6.390121
                |
treatment#month |
         TAU#0  |          0  (base)
         TAU#2  |          0  (base)
         TAU#3  |          0  (base)
         TAU#5  |          0  (base)
         TAU#8  |          0  (base)
       BtheB#0  |          0  (base)
       BtheB#2  |  -7.108148   2.924213    -2.43   0.015     -12.8395   -1.376796
       BtheB#3  |  -5.386667   2.924213    -1.84   0.065    -11.11802    .3446852
       BtheB#5  |  -4.318519   2.924213    -1.48   0.140    -10.04987    1.412833
       BtheB#8  |  -2.628148   2.924213    -0.90   0.369      -8.3595    3.103204
                |
          _cons |    10.7876   4.579099     2.36   0.018     1.812733    19.76247
---------------------------------------------------------------------------------

TAU群に比べたBtheB群のベースラインからの変化:meanCFB

contrast {treatment#month    1 -0.25 -0.25 -0.25 -0.25  /*baseTAU 2mTAU 3mTAU 5mTAU 8mTAU*/  ///
                                             -1  0.25  0.25  0.25  0.25} /*baseBtB 2mBtB 3mBtB 5mBtB 8mBtB*/


---------------------------------------------------
                |         df        chi2     P>chi2
----------------+----------------------------------
bdi             |
treatment#month |          1        3.31     0.0688
---------------------------------------------------

-----------------------------------------------------------------
                |   Contrast   Std. Err.     [95% Conf. Interval]
----------------+------------------------------------------------
bdi             |
treatment#month |
       (1) (1)  |   -4.86037   2.670517     -10.09449    .3737474
-----------------------------------------------------------------


主要評価項目を「治療期間全体での平均的なBDIの増減」とした場合,
BtheB群だけでみると,

・BtheB:mean

contrast {treatment#month   0    0    0    0    0  ///
                                           -1 0.25 0.25 0.25 0.25}


-----------------------------------------------------------------
                |   Contrast   Std. Err.     [95% Conf. Interval]
----------------+------------------------------------------------
bdi             |
treatment#month |
       (1) (1)  |  -12.12037   1.851671     -15.74958   -8.491163
-----------------------------------------------------------------


本に出てきたmean CFBのように設定すると,

contrast {treatment#month -1 0.25 0.25 0.25 0.25}

ですが,この結果は

-----------------------------------------------------------------
                |   Contrast   Std. Err.     [95% Conf. Interval]
----------------+------------------------------------------------
bdi             |
treatment#month |
       (1) (1)  |      -7.26   1.924313     -11.03158   -3.488417
-----------------------------------------------------------------


となって,SASの結果のpostの母数効果の値と一致するみたいです(p79参照)。
つまり,治療群に関係なく,ベースラインと比べた後の時点全体の効果。

おまけでp79 Model Vの再現

mixed bdi drug length treatment##post || id:post , reml  cov(unstruct) allbase


--------------------------------------------------------------------------------
           bdi |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+----------------------------------------------------------------
          drug |   1.123413   2.027048     0.55   0.579    -2.849528    5.096354
        length |   7.225293   1.988499     3.63   0.000     3.327906    11.12268
               |
     treatment |
          TAU  |          0  (base)
        BtheB  |  -1.816098   2.534022    -0.72   0.474    -6.782689    3.150493
               |
          post |
            0  |          0  (base)
            1  |      -7.26   1.924314    -3.77   0.000    -11.03159   -3.488413
               |
treatment#post |
        TAU#0  |          0  (base)
        TAU#1  |          0  (base)
      BtheB#0  |          0  (base)
      BtheB#1  |   -4.86037    2.67052    -1.82   0.069    -10.09449     .373752
               |
         _cons |   10.78761     4.5791     2.36   0.018     1.812743    19.76248
--------------------------------------------------------------------------------
0

複数変数の層別平均値を変数として追加し各層の最初の1行だけ表示

以前も以下の記事で紹介した

グループごとの平均値や割合を新たな変数として作成

について,さらに複数変数で定義されるグループ(層)の例
(例えば,学校AのBクラスなど)と,クラス内で全員が同じ値になる
集団平均値などの値を1つだけ取り出して扱いたいときなど(グラフを
描く時とか)に参考になるメモを書いておきます。

autoデータセットで例を示します。
まず,クロス表を出してみて,セルの人数を確認してみます。

sysuse auto
tab rep78 foreign


        Repair |
        Record |       Car type
          1978 |  Domestic    Foreign |     Total
    -----------+----------------------+----------
             1 |         2          0 |         2
             2 |         8          0 |         8
             3 |        27          3 |        30
             4 |         9          9 |        18
             5 |         2          9 |        11
    -----------+----------------------+----------
         Total |        48         21 |        69


rep78の1,2のforeignは0名なので,単純にするために今回はrep78の1,2の者を
削除しておきます。

drop if rep78==1 | rep78==2

rep78×foreignの層別でmpgの平均値を計算し,
新たな変数reppgを追加するには次のegenコマンドを使います。

egen reppg=mean(mpg), by(rep78 foreign)

結果を表示してみます。

sort rep78 foreign  //並び替え
list rep78 foreign mpg reppg


         | rep78    foreign   mpg      reppg |
      1. |     3   Domestic    17         19 |
      2. |     3   Domestic    20         19 |
                     (略)
     26. |     3   Domestic    19         19 |
     27. |     3   Domestic    14         19 |
     28. |     3    Foreign    23   23.33333 |
     29. |     3    Foreign    21   23.33333 |
     30. |     3    Foreign    26   23.33333 |
         |-----------------------------------|
     31. |     4   Domestic    14   18.44444 |
     32. |     4   Domestic    18   18.44444 |
     33. |     4   Domestic    22   18.44444 |
     34. |     4   Domestic    15   18.44444 |
     35. |     4   Domestic    14   18.44444 |
         |-----------------------------------|
     36. |     4   Domestic    28   18.44444 |
     37. |     4   Domestic    18   18.44444 |
     38. |     4   Domestic    16   18.44444 |
     39. |     4   Domestic    21   18.44444 |
     40. |     4    Foreign    23   24.88889 |
         |-----------------------------------|
     41. |     4    Foreign    23   24.88889 |
     42. |     4    Foreign    21   24.88889 |
     43. |     4    Foreign    25   24.88889 |
     44. |     4    Foreign    25   24.88889 |
     45. |     4    Foreign    24   24.88889 |
         |-----------------------------------|
     46. |     4    Foreign    28   24.88889 |
     47. |     4    Foreign    30   24.88889 |
     48. |     4    Foreign    25   24.88889 |
     49. |     5   Domestic    30         32 |
     50. |     5   Domestic    34         32 |
         |-----------------------------------|
     51. |     5    Foreign    18   26.33333 |
     52. |     5    Foreign    18   26.33333 |
                     (略)


このように,層別の平均値reppgが変数として追加されています。
ただこれではreppgの情報が重複していて,一覧で把握する際には結果が
多すぎて面倒です。なので,各層の最初の1行だけに表示したいと思います。

それには,以下のegenコマンドのtag関数で,各層の最初の1行を特定する変数
(ここではpickone)を作成すると便利です。

egen pickone = tag(rep78 foreign)

これで各層の最初の1行に1が割り当てられ,その他は0とする2値変数が
追加されたので,pickoneという変数の値が1のケースだけ表示すれば
よいわけです。

list rep78 foreign mpg reppg if pickone==1

         | rep78    foreign   mpg      reppg |
         |-----------------------------------|
      1. |     3   Domestic    17         19 |
     28. |     3    Foreign    23   23.33333 |
     31. |     4   Domestic    14   18.44444 |
     40. |     4    Foreign    23   24.88889 |
     49. |     5   Domestic    30         32 |
         |-----------------------------------|
     51. |     5    Foreign    18   26.33333 |
0

グラフの白黒化と凡例の位置の変更

autoデータを使って,車の値段をアウトカム,生産国(国産 or 外国),
重量,およびそれらの交互作用を説明変数とした線形回帰を実行します。

sysuse auto
reg price i.foreign##c.weight, nohead

--------------------------------------------------------------------------
           price|     Coef.  Std. Err.    t   P>|t|   [95% Conf. Interval]
----------------+---------------------------------------------------------
         foreign|
        Foreign | -2171.597  2829.409  -0.77  0.445  -7814.676    3471.482
          weight|  2.994814  .4163132   7.19  0.000   2.164503    3.825124
                |
foreign#c.weight|
        Foreign |  2.367227  1.121973   2.11  0.038    .129522    4.604931
                |
           _cons| -3861.719  1410.404  -2.74  0.008  -6674.681   -1048.757
--------------------------------------------------------------------------


次に,重量に応じた生産国ごとの値段の推定値(調整済み平均)を出します。
このためにはまず,連続量の重量をどのように区分するのか検討が必要です。

sum weight

 
   Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
     weight |        74    3019.459    777.1936       1760       4840



重量の最小値は1760,最大値が4840ということが分かります。
なので,とりあえず500ごとに区分してみようと思います。

margins foreign, at(weight=(1760(500)4840))

Adjusted predictions                              Number of obs   =    74
Model VCE    : OLS

Expression   : Linear prediction, predict()

1._at        : weight          =        1760

2._at        : weight          =        2260

3._at        : weight          =        2760

4._at        : weight          =        3260

5._at        : weight          =        3760

6._at        : weight          =        4260

7._at        : weight          =        4760

-----------------------------------------------------------------------
            |          Delta-method
            |   Margin   Std. Err.     t   P>|t|   [95% Conf. Interval]
------------+----------------------------------------------------------
_at#foreign |
1#Domestic  | 1409.153    708.814    1.99  0.051  -4.532181    2822.838
 1#Foreign  | 3403.875   727.8271    4.68  0.000    1952.27     4855.48
2#Domestic  |  2906.56   525.2356    5.53  0.000    1859.01    3954.109
 2#Foreign  | 6084.895   444.5963   13.69  0.000   5198.176    6971.614
3#Domestic  | 4403.966   368.7627   11.94  0.000   3668.492     5139.44
 3#Foreign  | 8765.915   639.0249   13.72  0.000    7491.42    10040.41
4#Domestic  | 5901.373   287.6764   20.51  0.000   5327.621    6475.126
 4#Foreign  | 11446.94   1077.865   10.62  0.000   9297.201    13596.67
5#Domestic  |  7398.78   340.8634   21.71  0.000   6718.949     8078.61
 5#Foreign  | 14127.96   1567.797    9.01  0.000   11001.08    17254.83
6#Domestic  | 8896.187   486.0826   18.30  0.000   7926.726    9865.648
 6#Foreign  | 16808.98   2072.905    8.11  0.000    12674.7    20943.25
7#Domestic  | 10393.59   665.5998   15.62  0.000   9066.097    11721.09
 7#Foreign  |    19490   2584.306    7.54  0.000   14335.76    24644.23
-----------------------------------------------------------------------

marginsを実行した直後に次のコマンドを打つと,グラフが出力されます。

marginsplot




論文に報告する時など,カラーではなく,白黒で分かるバージョンが
欲しかったりします。それはschemeオプションで指定すれば可能です。

marginsplot,   scheme(s1mono)



他の種類は

help scheme

で確認できます。


グラフの下部にある凡例を,グラフ領域の中に持っていきたい場合
legend()オプションをつけ,さらに中にring(0)を,そしてpos()で
位置を調整します。

marginsplot,   scheme(s1mono) legend(ring(0) pos(10))



pos()の中の数字は,時計の時刻の位置を表していて,
pos(10)は10時の位置,つまり左上の部分ということになります。

次に,横に広がっている凡例を,縦に表示するように
するには,col(1)を追加します。

marginsplot,   scheme(s1mono) legend(ring(0) pos(10) col(1))

0

ヒストグラムの基本2

グラフのタイトルなどを消してすっきりさせる方法を
メモします。

sysuse auto

まずは普通にmpgのヒストグラムを描きます

hist mpg,freq




縦軸と横軸のタイトルを削除します。" "の中に代わりに入れたい
タイトルを入れたら,それになります。

hist mpg,freq xtitle("") ytitle("")




次は層別ヒストグラム

hist mpg,freq by( foreign )




左下のGraphs by Car typeを消したい場合は,

hist mpg,freq by( foreign, note("") )




グラフの層をあらわすDomestic, Foreignの部分を消したい場合は,

hist mpg,freq by( foreign )  subtitle("")

0

ポリコリック相関行列を用いた探索的因子分析

解析にはStata13.1を使います。
次のリンク先を参考にしました。

Stata FAQ How can I perform a factor analysis with categorical (or categorical and continuous) variables?

データセットはRのパッケージpsychのbfiデータを使います。
bfiデータについてはこちら参照
単純にするために,変数A1~A5のみで解析します。
まずは欠損値の確認です。

misstable sum A1-A5


                                                     Obs<.
                                      +------------------------------
         |                            | Unique
Variable | Obs=.     Obs>.     Obs<.  | values        Min         Max
---------+----------------------------+------------------------------
      A1 |    16               2,784  |      6          1           6
      A2 |    27               2,773  |      6          1           6
      A3 |    26               2,774  |      6          1           6
      A4 |    19               2,781  |      6          1           6
      A5 |    16               2,784  |      6          1           6
---------------------------------------------------------------------

Obs=.の列を見ると,1変数につき16~27名で欠損値がみられることが
わかります。次にパターンをみます。

misstable patterns A1-A5, freq


      Missing-value patterns
        (1 means complete)

              |   Pattern
    Frequency |  1  2  3  4    5
  ------------+------------------
        2,709 |  1  1  1  1    1
              |
           22 |  1  1  1  1    0
           20 |  1  1  1  0    1
           15 |  1  0  1  1    1
           12 |  0  1  1  1    1
           12 |  1  1  0  1    1
            3 |  1  1  0  0    0
            2 |  0  1  0  1    1
            1 |  0  0  1  1    1
            1 |  0  1  1  0    1
            1 |  1  1  0  0    1
            1 |  1  1  0  1    0
            1 |  1  1  1  0    0
  ------------+------------------
        2,800 |

  Variables are  (1) A1  (2) A5  (3) A4  (4) A3  (5) A2

表の中身の1行目にあるように,変数のすべてで1であるもの,
つまり欠損値のないケースは2,709名だとわかります。

次にピアソンの相関行列をだします。

cor A1-A5


(obs=2709)
             |       A1       A2       A3       A4       A5
-------------+---------------------------------------------
          A1 |   1.0000
          A2 |  -0.3416   1.0000
          A3 |  -0.2683   0.4868   1.0000
          A4 |  -0.1484   0.3352   0.3622   1.0000
          A5 |  -0.1827   0.3878   0.5052   0.3067   1.0000


次はポリコリック相関です。そのままでは出ないので,
ユーザー作成コマンドのpolychoricをインストールする必要があります。
ネットにつながっている環境で

findit polychoric

と打ち,インストールしたら次のようにかきます。

polychoric A1-A5

Polychoric correlation matrix

            A1          A2          A3          A4          A5
A1           1
A2  -.40888172           1
A3  -.32566935    .5575669           1
A4  -.17597901   .38965116   .41069263           1
A5  -.22843963   .44653766   .57349574   .35401781           1

上記コマンドの実行で,様々な情報が生成されていて,
確認するには次のコマンドを打ちます。

return list

scalars:
               r(pLR0) =  6.15755687868e-62
                r(LR0) =  275.8067214420807
                r(pX2) =  .0003008227338877
               r(dfX2) =  24
                 r(X2) =  55.1282052450148
                r(pG2) =  .0012040446647921
               r(dfG2) =  24
                 r(G2) =  50.55156075192645
             r(se_rho) =  .0201740996318403
                r(rho) =  .3540178066269911
                  r(N) =  2709
              r(sum_w) =  2709

macros:
               r(type) : "polychoric"

matrices:
                  r(R) :  5 x 5

これらの情報は一時的なので次の解析を行うと消えてしまいます。
解析に用いた人数など,個別の値だけを見たいときには,
以下のようにします。

display r(sum_w)

2709

※r(N)とr(sum_w)の違いはヘルプを見ても詳しく書いてなかったので,
不明です


そして,この情報を用いた次の3行がポリコリック相関行列を
用いた探索的因子分析です。3行をいっぺんに実行します。

local N = r(sum_w)
matrix r = r(R)
factormat r, n(
`N') factors(1)

 Factor analysis/correlation                Number of obs    =   2709
   Method: principal factors                Retained factors =      1
   Rotation: (unrotated)                    Number of params =      5

   --------------------------------------------------------------------
        Factor  |   Eigenvalue  Difference   Proportion   Cumulative
   -------------+------------------------------------------------------
       Factor1  |      1.95127     1.85836       1.1904       1.1904
       Factor2  |      0.09291     0.14472       0.0567       1.2470
       Factor3  |     -0.05181     0.12276      -0.0316       1.2154
       Factor4  |     -0.17457     0.00401      -0.1065       1.1089
       Factor5  |     -0.17858           .      -0.1089       1.0000
   -----------------------------------------------------------------------
 LR test: independent vs. saturated:  chi2(10) = 3390.61 Prob>chi2 = 0.0000

Factor loadings (pattern matrix) and unique variances

   ---------------------------------------
       Variable |  Factor1 |   Uniqueness
   -------------+----------+--------------
             A1 |  -0.4376 |      0.8085 
             A2 |   0.7078 |      0.4990 
             A3 |   0.7558 |      0.4287 
             A4 |   0.5154 |      0.7343 
             A5 |   0.6494 |      0.5782 
   ---------------------------------------

コマンドの解説をすると,まずlocalマクロ変数でポリコリック相関行列算出により出てきた
人数r(sum_w)をNと定義しています。次にポリコリック相関行列r(R)をrと定義します。
そして,3行目が相関行列を用いた探索的因子分析です。
因子数は1に指定しています。factormatコマンドの後に相関行列rを置き,オプションで
相関行列を出すのに用いた人数Nをいれています。デフォルトでは主因子法*です。

次にピアソン相関行列を用いた探索的因子分析をしてみます。
先ほど出したので,相関行列の表示は省略することにし,先頭にquiをつけます。

qui cor A1-A5
return list


scalars:
                  r(N) =  2709
                r(rho) =  -.3416241810153322

matrices:
                  r(C) :  5 x 5

今度は,格納されている結果の名前が少し違っています。

local N = r(N)
matrix r = r(C)
factormat r, n(`N') factors(1)

Factor analysis/correlation                   Number of obs=     2709
Method: principal factors                 Retained factors =        1
Rotation: (unrotated)                     Number of params =        5

---------------------------------------------------------------------
Factor       Eigenvalue   Difference   Proportion   Cumulative
-------------+-------------------------------------------------------
Factor1        1.66236      1.59574       1.2699       1.2699
Factor2        0.06662      0.13026       0.0509       1.3208
Factor3       -0.06364      0.10961      -0.0486       1.2722
Factor4       -0.17325      0.00978      -0.1323       1.1398
Factor5       -0.18303            .      -0.1398       1.0000
----------------------------------------------------------------------
LR test: independent vs. saturated:  chi2(10) = 2531.30 Prob>chi2 = 0.0000

Factor loadings (pattern matrix) and unique variances

---------------------------------------
Variable   Factor1    Uniqueness
-------------+----------+--------------
A1       -0.3868       0.8504 
A2       0.6509       0.5763 
A3       0.7029       0.5059 
A4       0.4813       0.7684 
A5       0.6027       0.6367 
---------------------------------------

大枠は同じですが,数値はかなり違ったようになってきます。
ちなみに,ピアソン相関行列で行った因子分析は,下記のように
基本の因子分析を行った場合と同じ結果になります。

factor A1-A5, factors(1)



*推定法について,清水裕士先生(@simizu706)の↓のツイート

ポリコリック相関行列をそのまま因子分析に適用するのは,あまりよくない。
標準誤差を重みづけた,重みづき最小二乗法(WLS)を用いて推定したほうがいい。

が気になるので,この辺についてもその内Stataでのやり方を調べてみるつもりです。
SEMではできるぽいので。

本記事の内容について,清水先生のブログ記事↓も参考になります

カテゴリカルデータの相関係数
0