Stataメモ
Stataメモ目次
投稿日時 : 2016/05/05
土屋 政雄
土屋 政雄
参照に便利なように，目次を作りました
researchmapのバージョンアップによるリンク切れを修正 2020223
【データハンドリング】
- データハンドリングのよく使うコマンド集（自分用）
- SPSSのデータセットをStataで開く
- 変数ラベルと値のラベル
- 連続量の変数からカテゴリ変数を作成（カットオフ分割，ダミー変数）
- 文字型の日付データを日付変数に変換して計算できるようにする
- 2つの時刻の時間と分が別々の列にある変数群から経過時間変数を作成
- たくさんの変数に同じ処理を繰り返したいとき
- 連続して変わる数値の繰り返し処理を入れ子で実行する（double loop）
- 変数の複数の値のいずれかに該当したら実行するという条件の書き方
- 変数を足し合わせた得点で欠損があるときの平均値の計算
- グループごとの平均値や割合を新たな変数として作成
- 複数変数の層別平均値を変数として追加し各層の最初の1行だけ表示
- 層（グループ）別にidと合計人数を出し変数として追加（_nと_N）
- データセットをwideからlongへ変換する
- 追跡データの欠損について前の時点の値を代入する（LOCF）
【解析】
- 文字型変数の値が長すぎてtabulateの結果で文字がとぎれるのを解決
- 層別に平均値と標準偏差（SD）と人数だけを算出
- 層別に複数変数の平均値を縦一列に算出
- summarizeで要約統計量を出すときに出力に変数ラベルを表示させる
- データの欠損を数える
- 平均値，標準偏差，人数の値だけからt検定を行う
- 割合の95%信頼区間を算出する
- t検定と効果量（エフェクトサイズ）
- 相関係数の95%信頼区間
- 相関行列を出す際のpwcorrとcorrelateの違い
- 点双列相関係数の出し方
- 再検査（再テスト）信頼性のための級内相関係数(ICC)(2):Stataで算出
- 分散分析と線形回帰のオムニバス効果量
- 線形回帰のメモ （いわゆる単回帰分析，重回帰分析）
- 線型回帰分析におけるカテゴリ×連続の交互作用
- 階層的重回帰
- ロジスティック回帰のメモ
- 複数の回帰分析の結果をtableで報告する形式に出力する
- 介入研究の解析での線形混合モデル（1） 修正
- 介入研究の解析での線形混合モデル（3）介入効果はどこを見るか 修正
- ROC曲線を描く
- 解析結果をExcelなどへ出力する（2要因のクロス表など）
- ポリコリック相関行列を用いた探索的因子分析
- 検証的（確証的／確認的）因子分析：confirmatory factor analysis(CFA)
- 検証的（確証的／確認的）因子分析 (2)因子間相関
【図・グラフ】
0
複数の回帰分析の結果をtableで報告する形式に出力する
投稿日時 : 2016/04/10
土屋 政雄
土屋 政雄
回帰分析で説明変数が1つだけの2変量モデルから，調整する変数
を段階的に加えていき，複数のモデルを並べてtableで報告することが
よく行われます。普通にやったら手動で編集して作表しなければ
なりませんが，ユーザー作成コマンドのestoutを使うととても
楽にできます。
まずパッケージestoutをインストールします。まだ入ってない場合，
ネットにつながっている状態で次のように打ちます。
(1)回帰係数のみを出力
まずはアウトカムをforeign，説明変数をmpgだけにしたロジスティック回帰
ここで，次のコマンドを実行すると，回帰係数だけが出てきます。
ただし，オッズ比ではありません。
念のため，オッズ比変換前の回帰係数も出してみて確認
(2)複数の回帰分析の結果を出力
さて，それでは，1つずつ説明変数を加えて行った，計3つのモデルを
作成します。出力結果を出さないように，quietlyを追加しています。
これで3つのモデルが格納されたので確認してみます。
3つのモデルすべてを出力します。*はすべてのモデルを指定する記号です。
一列目が説明変数で，モデルごとに列が変わり，説明変数が増えていく
様が分かります。
(3)出力結果を整える
出したい形式は，コマンドが複雑なので，ちょっとずつ
オプションを追加して説明していきます。
・回帰係数と95%信頼区間を表示し，切片を削除
それぞれの回帰係数の下の行に信頼区間が追加されています
・オッズ比とオッズ比の95%信頼区間
・小数点第1位まで表示し，切片を削除
・95%信頼区間を横に配置
bからciの部分を" "で囲みます
・95%信頼区間を[ ]で閉じる
おそらくこれが一番有用な形式と思います。ciの所にparを追加します
・95%信頼区間を( )で閉じる
parの所に( , )を追加します。
を段階的に加えていき，複数のモデルを並べてtableで報告することが
よく行われます。普通にやったら手動で編集して作表しなければ
なりませんが，ユーザー作成コマンドのestoutを使うととても
楽にできます。
まずパッケージestoutをインストールします。まだ入ってない場合，
ネットにつながっている状態で次のように打ちます。
ssc install estout
autoデータセットを使います。
rep78*foreignのクロス表では，0名のセルが発生するので，
単純にするためにrep78の3未満をdropしておきます。
sysuse auto
drop if rep78<3
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
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）
投稿日時 : 2016/02/27
土屋 政雄
土屋 政雄
本記事の解説では，Stata 13.1を使います。
ここで解説するRCTなどの介入研究で使われる解析法は，
線形混合モデル（linear mixed model），マルチレベル(multilevel)，
階層線型モデル(hierarchical linear model)など様々な呼び方が
ありますが，ここでは介入研究の文脈でよく用いられる線形混合モデル
の呼び方にします（「線型」という表記もあります）。
従来は，アウトカムが連続量であるRCTなどの3時点以上の全体の解析は
反復測定分散分析（repeated ANOVA）で行われていました。しかし近年では，
これが線形混合モデルになることが多いので，こうしたタイプの
解析を行う研究者にとって必須の知識になりつつあると思います。
特に，追跡時に必ずと言っていいほど発生する欠損について，
データがそろっている対象者だけに限定するとか，前の時点の値を代入する
（LOCF）とかをしなくても，問題なく解析できる点が強みです。
LOCFについては↓参照
いろいろ勉強しないとわからないことが多いのですが，まず解析
して体験することでわかることも多いと思いますので，
これまで解説に使ってきたBtheBデータを使って試してみます。
ちなみにBtheBデータセットは，実際の研究のデータで，↓のように
論文も出てるので，解釈の参考にもなります。
BtheBデータの使い方は下記のページ↓参照
BtheBデータにある欠損値の状況は下記のページ↓で解説してます
サンプルデータで現実的な欠損値が含まれているものは少ない気がする
ので，これがこのデータセットの強みです。
線型混合モデルで縦断データを解析する際の準備は，下記のページの作業↓
が必要になります
さて，これで準備が整いましたので，いよいよ解析に進みます。
基本的な解析自体は非常に簡単で，抑うつ症状（bdi）について
群(treatment)×時間(month)の交互作用の検証という主目的は以下の
1行だけです。以下はランダム切片だけですが，Stataだとこんなシンプルに
書けるのが素敵です。
##の記号で説明変数をつなぐだけで，
それぞれの主効果も自動で含めてくれます。Stataではバージョン12まではxtmixedという
コマンドでしたが，13からmixedになっており，徐々に進化している感があります。
詳しくはhelp mixedと打つといろいろわかります。
出力で最も関心のある交互作用を確認するには↓のtreatment#monthの所を見ます。
※2016/01/31 修正 オムニバスの交互作用は最近ではあまり使われないです。詳しくは
介入研究の解析での線形混合モデル（3）介入効果はどこを見るか 参照
----------------------------------------------------------------------------------------------------
以下は記事修正前に書いていた，オムニバス交互作用の検定をしたことを前提に示した
結果となります。
デフォルトだとmonthがカテゴリ変数として認識されていました。
これだけでは，論文に記載することの多い交互作用項全体のp値はわかりません。
これを出すためには，以下のコマンドを実行します。
これによって，群×時間の交互作用はp=0.3266で有意ではないという
結果が分かりました。
このカイ二乗値をF値にしたいときは，↓を参考に
カイ二乗値を自由度で割る計算をするだけです。すなわち
というように上のカイ二乗値を自由度で割るコマンドを書くと，
という結果が出てきます。これがF値になります。
ちなみに，Stata公式がカイ二乗値とF値の関係について解説しています
入門知識はこちらにまとめていますのでご参考まで
関連記事：
このトピックスは適宜修正・追加予定です。
【履歴】
2014/06/02 記事作成
2014/07/21 一部削除→交互作用の部分
2016/01/31 修正。出力のポイントについて説明と参照リンク追加
2016/02/27 前回修正した部分の表現を修正
ここで解説するRCTなどの介入研究で使われる解析法は，
線形混合モデル（linear mixed model），マルチレベル(multilevel)，
階層線型モデル(hierarchical linear model)など様々な呼び方が
ありますが，ここでは介入研究の文脈でよく用いられる線形混合モデル
の呼び方にします（「線型」という表記もあります）。
従来は，アウトカムが連続量であるRCTなどの3時点以上の全体の解析は
反復測定分散分析（repeated ANOVA）で行われていました。しかし近年では，
これが線形混合モデルになることが多いので，こうしたタイプの
解析を行う研究者にとって必須の知識になりつつあると思います。
特に，追跡時に必ずと言っていいほど発生する欠損について，
データがそろっている対象者だけに限定するとか，前の時点の値を代入する
（LOCF）とかをしなくても，問題なく解析できる点が強みです。
LOCFについては↓参照
追跡データの欠損について前の時点の値を代入する（LOCF）
いろいろ勉強しないとわからないことが多いのですが，まず解析
して体験することでわかることも多いと思いますので，
これまで解説に使ってきたBtheBデータを使って試してみます。
ちなみにBtheBデータセットは，実際の研究のデータで，↓のように
論文も出てるので，解釈の参考にもなります。
一般診療の場における不安と抑うつへのコンピュータ化された，
対話式の，マルチメディア認知行動療法プログラム
Proudfoot et al., Psychol Med. 2003 ;33(2):217-27.
プライマリケアにおける不安と抑うつへのコンピューター認知行動療法の
臨床的効果：無作為割付比較試験
Proudfoot et al. Br J Psychiatry. 2004 ;185:46-54.
対話式の，マルチメディア認知行動療法プログラム
Proudfoot et al., Psychol Med. 2003 ;33(2):217-27.
プライマリケアにおける不安と抑うつへのコンピューター認知行動療法の
臨床的効果：無作為割付比較試験
Proudfoot et al. Br J Psychiatry. 2004 ;185:46-54.
BtheBデータの使い方は下記のページ↓参照
EZRでパッケージ中のRデータセットをStataデータセットに変換する
BtheBデータにある欠損値の状況は下記のページ↓で解説してます
サンプルデータで現実的な欠損値が含まれているものは少ない気がする
ので，これがこのデータセットの強みです。
データの欠損を数える
線型混合モデルで縦断データを解析する際の準備は，下記のページの作業↓
が必要になります
データセットをwideからlongへ変換する
さて，これで準備が整いましたので，いよいよ解析に進みます。
基本的な解析自体は非常に簡単で，抑うつ症状（bdi）について
群(treatment)×時間(month)の交互作用の検証という主目的は以下の
1行だけです。以下はランダム切片だけですが，Stataだとこんなシンプルに
書けるのが素敵です。
mixed bdi treatment##month || id:
##の記号で説明変数をつなぐだけで，
それぞれの主効果も自動で含めてくれます。Stataではバージョン12まではxtmixedという
コマンドでしたが，13からmixedになっており，徐々に進化している感があります。
詳しくはhelp mixedと打つといろいろわかります。
出力で
※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
割合と信頼区間のグラフを描く
投稿日時 : 2016/02/27
土屋 政雄
土屋 政雄
割合の95%信頼区間を算出する
にて，割合と信頼区間の算出を紹介しました。
これらをStataでをグラフにするには，
割合と信頼区間それぞれが変数となった新たなデータセットが必要です。
以下，バージョンはStata14.1を用いました。
上記の解析で得られた結果出力をexcel等で加工して取り出して
作ることもできますが，手間がかかります。
そこで，何とかStata上だけで処理する方法を考えてみました。
1. 割合と信頼区間の情報のみを抽出
例えば，rep78=3だけ取り出してみてみます。
この解析が行われると，いくつかの情報がメモリに格納されるようです。
と打つと出てきます。
この内，割合と信頼区間の情報である
さえ選んで出力できれば便利です。以下のようにしてできます。
出力したい値をコンマ（,）で区切って並べると，スペースで
区切って出力されます
あとは，これをrep78の値ごとに繰り返せば
便利です。そんな時はloopを使います。くわしくは
→たくさんの変数に同じ処理を繰り返したいとき
rep78が3-5まであるので，それぞれごとに
foreignの割合と信頼区間を算出し，算出結果の
割合，95%信頼区間の下限，上限を，rep78の値ごと（i）
にdisplayで選んで出力。なお，ciの出力はquietlyで抑制
出力は以下の通り
あとはこれをStataデータセットに読み込めれば，
グラフ化に便利です。
2. 新たなデータセット作成
ここでいったんデータセットをclearします。
1.で取り出した情報を結果画面からコピーして
以下のようにエディタに書きます。
inputの行は変数の名前です。そして最後にendを置きます。
次に，100%換算にするために，割合と信頼区間
それぞれに100を掛けます。
そして，新たなデータセットを確認します。
これで準備完了です。
3. 割合と信頼区間のグラフ作成（縦書き）
ここまでできたらあとは簡単です。
rspikeで信頼区間，scatterで割合をそれぞれ描画します。
見栄えが悪いので，オプションを色々追加してグラフを表示します。
msymbolで点推定値を四角形に指定し，msizeで大きさを変更，
xlabelで，rep78が3～5まで1ずつ表示されるようにし，
legendで凡例の位置を調整，ラベルも変更，plotregionで左右の空きを増やしました。
rspikeをrcapに変えると，信頼区間の端に棒がつきます。
4. 割合と信頼区間のグラフ作成（横書き）
横向きに描きたいときは，rcapの方にhorizontalオプションを付け，
scatterの変数の順序を逆にします。
にて，割合と信頼区間の算出を紹介しました。
これらを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)
}
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
/*↓に貼り付け*/
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
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(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))
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%信頼区間を算出する
投稿日時 : 2016/02/24
土屋 政雄
土屋 政雄
stata14.1を使います。使用するコマンドci proportionsは
バージョン14.1から追加されたみたいです。
以前はciでbinomialオプションで扱っていたようです。
autoデータを使います。
複数変数の層別平均値を変数として追加し各層の最初の1行だけ表示
の記事でみたように，rep78とforeignのクロス表では0セルが
発生するので，rep78の1と2は削除しておきます。
1.割合の95%信頼区間
まず変数foreignを使います。
0:Domestic, 1:Foreignという値なので，
1:Foreignの割合を算出するということになります。
まずは確認で度数分布
rep78の1と2をdropしているので人数は減っています。
次に割合と信頼区間
proportionが0.34なので，34%です。信頼区間は0.229, 0.473
defaultの信頼区間はexact（Clopper-Pearson）で，
オプションでwald, wilson, agresti, jeffreysに対応しています
2.層別に割合の95%信頼区間
rep78の値ごとに出したい場合は
最後にrep78が欠損の場合も出力されています
バージョン14.1から追加されたみたいです。
以前はciでbinomialオプションで扱っていたようです。
autoデータを使います。
複数変数の層別平均値を変数として追加し各層の最初の1行だけ表示
の記事でみたように，rep78とforeignのクロス表では0セルが
発生するので，rep78の1と2は削除しておきます。
sysuse auto
drop if rep78==1 | rep78==2
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）介入効果はどこを見るか
投稿日時 : 2016/01/31
土屋 政雄
土屋 政雄
線形混合モデルを使用した際の、治療効果や介入効果の求め方について、
以下の本，
経時的繰り返し測定デザイン : ─治療効果を評価する混合効果モデルとその周辺─ (医学統計学シリーズ)
丹後 俊郎
朝倉書店(2015/09/19)
値段:￥ 4,860
が非常に詳細に解説していて，とても理解が進んでよかったので、
本の中ではSASで示されていたコード部分の例をStataで再現してみます。
Stataのバージョンは13.1です。
BtheBデータのStataへの読み込み方はたとえば↓の過去記事参照
読み込んだらまず，後にデータセット構造をwideからlongに変換する必要が
あるので，Stataのルール上，変数名を次のように変換します。
次にIDがないので作成します。
最初の5名のデータを確認します。
(1)アウトカムに１つでも欠損のある者を除いた完全ケースデータを作成（p74のModel II）
idの1,3,5番は欠損値があったので，データセットから除かれています。
次に群ごとの人数を確認します。
p68にある群ごとの人数と一致します。
(2)データセットをwideからlongに変換
(3)線形混合モデル（p74のModel II） ※2016/01/31 本の解析と合うように修正
本ではベースライン以外の測定時を1とした変数postを
作成しているので同様に作ります。
本と同様に，treatmentを連続変数として扱うため以下の
様に値を0,1に変換します。
次が線形混合モデル本体です。
allbaseオプションは，referenceカテゴリーが何か明示するためにつけてます。
treatmentは連続変数として扱うため，c.をつけています。
カテゴリ変数のreferenceは，Stataではデフォルトが最初のカテゴリなので，SAS
のように変更する必要はありません。
ランダム効果部分の出力
母数効果の部分の出力
この各交互作用項の推定値が，TAU群に比べたBtheB群のベースラインからの変化について
の各測定時点での群間差になり，つまり効果量（非標準化）ということになります。
つまり，大抵の論文でのメインで報告する値です。どの測定時を主要評価時点
とするかは事前に設定しておく必要があります。
TAU群に比べたBtheB群のベースラインからの変化（介入期間を通じた平均）：meanCFB
--------------------------------------------------------------------------------------------------------
以下はtreatmentがカテゴリ変数だった場合の結果（2016/01/31以前の記事で説明
していた内容を要約）
ランダム効果部分の出力
母数効果の部分の出力
TAU群に比べたBtheB群のベースラインからの変化：meanCFB
主要評価項目を「治療期間全体での平均的なBDIの増減」とした場合，
BtheB群だけでみると，
・BtheB:mean
本に出てきたmean CFBのように設定すると，
ですが，この結果は
となって，SASの結果のpostの母数効果の値と一致するみたいです（p79参照）。
つまり，治療群に関係なく，ベースラインと比べた後の時点全体の効果。
おまけでp79 Model Vの再現
以下の本，
経時的繰り返し測定デザイン : ─治療効果を評価する混合効果モデルとその周辺─ (医学統計学シリーズ)丹後 俊郎
朝倉書店(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)アウトカムに１つでも欠損のある者を除いた完全ケースデータを作成（p74のModel II）
keep if bdi0~=. & bdi2~=. & bdi3~=. & bdi5~=. & bdi8~=.
list in 1/5
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
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*/
-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}
-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行だけ表示
投稿日時 : 2015/11/25
土屋 政雄
土屋 政雄
以前も以下の記事で紹介した
グループごとの平均値や割合を新たな変数として作成
について，さらに複数変数で定義されるグループ（層）の例
（例えば，学校AのBクラスなど）と，クラス内で全員が同じ値になる
集団平均値などの値を１つだけ取り出して扱いたいときなど（グラフを
描く時とか）に参考になるメモを書いておきます。
autoデータセットで例を示します。
まず，クロス表を出してみて，セルの人数を確認してみます。
rep78の1,2のforeignは0名なので，単純にするために今回はrep78の1,2の者を
削除しておきます。
rep78×foreignの層別でmpgの平均値を計算し，
新たな変数reppgを追加するには次のegenコマンドを使います。
結果を表示してみます。
このように，層別の平均値reppgが変数として追加されています。
ただこれではreppgの情報が重複していて，一覧で把握する際には結果が
多すぎて面倒です。なので，各層の最初の1行だけに表示したいと思います。
それには，以下のegenコマンドのtag関数で，各層の最初の1行を特定する変数
（ここではpickone）を作成すると便利です。
これで各層の最初の1行に1が割り当てられ，その他は0とする2値変数が
追加されたので，pickoneという変数の値が1のケースだけ表示すれば
よいわけです。
グループごとの平均値や割合を新たな変数として作成
について，さらに複数変数で定義されるグループ（層）の例
（例えば，学校AのBクラスなど）と，クラス内で全員が同じ値になる
集団平均値などの値を１つだけ取り出して扱いたいときなど（グラフを
描く時とか）に参考になるメモを書いておきます。
autoデータセットで例を示します。
まず，クロス表を出してみて，セルの人数を確認してみます。
sysuse auto
tab rep78 foreign
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
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
グラフの白黒化と凡例の位置の変更
投稿日時 : 2015/09/20
土屋 政雄
土屋 政雄
autoデータを使って，車の値段をアウトカム，生産国（国産 or 外国），
重量，およびそれらの交互作用を説明変数とした線形回帰を実行します。
次に，重量に応じた生産国ごとの値段の推定値（調整済み平均）を出します。
このためにはまず，連続量の重量をどのように区分するのか検討が必要です。
重量の最小値は1760，最大値が4840ということが分かります。
なので，とりあえず500ごとに区分してみようと思います。
marginsを実行した直後に次のコマンドを打つと，グラフが出力されます。
論文に報告する時など，カラーではなく，白黒で分かるバージョンが
欲しかったりします。それはschemeオプションで指定すれば可能です。
他の種類は
で確認できます。
グラフの下部にある凡例を，グラフ領域の中に持っていきたい場合
leｇend()オプションをつけ，さらに中にring(0)を，そしてpos()で
位置を調整します。
pos()の中の数字は，時計の時刻の位置を表していて，
pos(10)は10時の位置，つまり左上の部分ということになります。
次に，横に広がっている凡例を，縦に表示するように
するには，col(1)を追加します。
重量，およびそれらの交互作用を説明変数とした線形回帰を実行します。
sysuse auto
reg price i.foreign##c.weight, nohead
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
で確認できます。
グラフの下部にある凡例を，グラフ領域の中に持っていきたい場合
leｇend()オプションをつけ，さらに中に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
投稿日時 : 2015/06/21
土屋 政雄
土屋 政雄
グラフのタイトルなどを消してすっきりさせる方法を
メモします。
まずは普通にmpgのヒストグラムを描きます
縦軸と横軸のタイトルを削除します。" "の中に代わりに入れたい
タイトルを入れたら，それになります。
次は層別ヒストグラム
左下のGraphs by Car typeを消したい場合は，
グラフの層をあらわすDomestic, Foreignの部分を消したい場合は，
メモします。
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
ポリコリック相関行列を用いた探索的因子分析
投稿日時 : 2015/06/14
土屋 政雄
土屋 政雄
解析には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)
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
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)
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