link管理
链接快照平台
- 输入网页链接，自动生成快照
- 标签化管理网页链接
相关文章推荐
|
时尚的熊猫 · Stata编程技巧(3): ...· 2 周前 · |
|
|
阳刚的火腿肠 · 如何在stata中将字符型数据库 | ...· 2 周前 · |
|
|
纯真的柠檬 · Stata数值型变成字符型_千锋教育· 2 周前 · |
|
|
一身肌肉的山楂 · Stata学习笔记|数据处理3-类型转化_学 ...· 2 周前 · |
|
|
豪爽的灯泡 · 《Stata软件操作》——经济学拔尖基地举办 ...· 2 周前 · |
|
|
越狱的刺猬 · 快穿主受文集百度云 - 百度· 8 月前 · |
|
|
爱跑步的电池 · 盘点86版《西游记》已经去世的17位演员，他 ...· 1 年前 · |
|
|
沉稳的茶壶 · random --- 生成伪随机数 — ...· 1 年前 · |
|
|
健身的咖啡豆 · Fortran程序中直接使用度（°）计算三角 ...· 1 年前 · |
|
|
不拘小节的香烟 · The process started ...· 1 年前 · |
【図・グラフ】
を段階的に加えていき，複数のモデルを並べて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本記事の解説では，Stata 13.1を使います。
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)
ここで解説する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 前回修正した部分の表現を修正 割合の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)
}
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
stata14.1を使います。使用するコマンドci proportionsは
バージョン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が欠損の場合も出力されています 線形混合モデルを使用した際の、治療効果や介入効果の求め方について、
以下の本，
経時的繰り返し測定デザイン : ─治療効果を評価する混合効果モデルとその周辺─ (医学統計学シリーズ)
丹後 俊郎
朝倉書店(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
--------------------------------------------------------------------------------
グループごとの平均値や割合を新たな変数として作成
について，さらに複数変数で定義されるグループ（層）の例
（例えば，学校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 |autoデータを使って，車の値段をアウトカム，生産国（国産 or 外国），
|-----------------------------------|
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 |
重量，およびそれらの交互作用を説明変数とした線形回帰を実行します。
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))
グラフのタイトルなどを消してすっきりさせる方法を
メモします。
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("")
解析には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)
※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ではできるぽいので。
本記事の内容について，清水先生のブログ記事↓も参考になります
カテゴリカルデータの相関係数
相関係数を出すときに使うpwcorr，およびcorrelate（略してcor）は，
どちらも指定した変数群の相関行列を出してくれますが，欠損値の有無で
結果が違います。例えば，Rのパッケージpsychのbfiデータの
変数A1-A5で試してみます。bfiデータについては
こちら
参照
pwcorr A1-A5
| A1 A2 A3 A4 A5
-------------+---------------------------------------------
A1 | 1.0000
A2 | -0.3402 1.0000
A3 | -0.2652 0.4851 1.0000
A4 | -0.1464 0.3351 0.3604 1.0000
A5 | -0.1814 0.3901 0.5041 0.3075 1.0000
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
この違いは，単純に計算に使われた人数の違いです。
corの方は，欠損のあるケースをリストワイズするからです。
そうすると(obs=2709)になるということが結果にも出ています。
したがって，pwcorrをリストワイズで出すオプションをつけると，
pwcorr A1-A5, list
| 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
となり，corの結果と一致します。
pwcorrのデフォルトではどのように解析がされているか確認するには，obs
オプションをつけます
pwcorr A1-A5, obs
| A1 A2 A3 A4 A5
-------------+---------------------------------------------
A1 | 1.0000
| 2784
|
A2 | -0.3402 1.0000
| 2757 2773
|
A3 | -0.2652 0.4851 1.0000
| 2759 2751 2774
|
A4 | -0.1464 0.3351 0.3604 1.0000
| 2767 2758 2759 2781
|
A5 | -0.1814 0.3901 0.5041 0.3075 1.0000
| 2769 2757 2758 2765 2784
相関係数の下にそれぞれのペアの人数が出てきました
Double loops
まずは単純な例。iの部分が1から3まで変わっていき，
計算結果がそれぞれ出力されます。
forvalues i=1/3{
display 1+`i'
}
display 1+`i'
}
2
3
4
このように1+1，1+2，1+3の計算結果が出力されました。
これを，式の部分もテキストとして表示したいときには，
localマクロを使います。式の部分はiがそのまま変わる
だけですが，答えの部分がiとは異なるためです。ここでは
答えの部分をsumと置いています。
forvalues i=1/3{
local sum=1+`i'
display "1+`i'=`sum'"
}
local sum=1+`i'
display "1+`i'=`sum'"
}
1+1=2
1+2=3
1+3=4
次は，入れ子構造のloopです。具体的には，式の1+の部分も数値を
増やしていきたいときの場合です。この方法は，double loopまたは
nested loopと呼ばれるみたいです。
forvalues i=1/3{
forvalues j=1/3{
local sum=`i'+`j'
display "`i'+`j'=`sum'"
}
}
forvalues j=1/3{
local sum=`i'+`j'
display "`i'+`j'=`sum'"
}
}
1+1=2
1+2=3
1+3=4
2+1=3
2+2=4
2+3=5
3+1=4
3+2=5
3+3=6
関連記事として， たくさんの変数に同じ処理を繰り返したいとき も参照
変更履歴
20150507 題名修正 先日まとめた
再検査（再テスト）信頼性のための級内相関係数(ICC)
について，今度はどうやってソフトで出すのかをまとめます。
前回と同じく，内容には自信がないので信用しないでください。
自分の把握している限り，Stata13の時点ではICCを出すのに2つの専用コマンドが
あります（線形混合モデルコマンドでも出せるみたいですが，ここでは扱いません）。
■ icc23
まずは， Paulsen et al. (2012) に書かれていた，従来Stataで使われてきた
コマンド"icc23"からみてみます。
これはユーザー作成コマンドなので，まずネットにつながっている環境で
インストールする必要があります。
findit icc23
以下，ヘルプより
help icc23
このコマンドは反復測定ANOVAに基づき，ICC[2,1], ICC[2,k], ICC[3,1],
and ICC[3,k]についてrandom effects modelsのICCを算出するものです。
データはまずreshapeで long形式にする 必要があり，その後以下のコマンド例に変数を
あてはめていく感じです。
icc23 score rater person_id
つまり，変数の順として， 尺度の得点 評定者 被評定者ID ，という並びです。
ただ，これでは評定者間のICC例なので，いまいちぴんときません。
いろいろ探した結果， 解析例 を見つけました。
icc23 pf time id, model(3)
つまり，time（測定ポイント）が評定者（rater）の位置にくるだけみたいです。
これは，下記の本の解析例のようです。
Data and Programming for Applied Categorical and Count Data Analysis
本を見たら，詳しく記述されているかもしれません。
ここではmodel(3)オプションをつけることで， ICC[3,1] または ICC[3,k]にしている
ようです。デフォルトではICC[2,1] または ICC[2,k]になります。
多くの場合，[ ,k]の方は使われないそうです。
こうした知識を入れた上で，icc23のヘルプを見直したら，ちゃんと書いてありました。
icc23 <dv> <classvar> <within_var>
classvarとは，個人内で繰り返される因子であり，
例として評定者，機器， 測定ポイント などがある
within_varとは，"被験者内"変数であり，例として，
評定される個人のことを指す
注：赤字は土屋が強調した
■ icc
一方，Stata13（12からかも）では，ICCの公式コマンドが追加されています。
Intraclass correlation coefficients
これも，ヘルプの解析例が評定者間なので，いまいちわかりません。
こちらもいろいろ探した結果，Chuck Huber氏のスライド「Psychometrics Using Stata」
に例がありました（ http://www.stata.com/meeting/sandiego12/abstracts/
からスライド，データやプログラムがDL可能）。以下がコマンド例です。
icc qu1_t id time
注意しないといけなそうなのは， idの位置がicc23と逆になる と思われる点です。
これらの位置をかえると結果が全然違うものになってしまいます。
解析の詳しい内容は公開されているdoファイルの中に書いてあるのですが，
要約すると，Huber_2012SanDiego_Pilot.dtaを読み込んで，reshapeで
long形式に直し，上記のコマンドを走らせるという手順です。
解釈としては，スライドで親切に赤字で示されている通り，
individualのICCの所が該当するということです。たぶんこれが(2,1)。
averageの方が(2,k)のような気がします。
↓
/*（20150429 追記）
Stataのマニュアルに詳しい記載がありました。
STATA BASE REFERENCE MANUAL RELEASE 13のp857
→題名でググればStata社のサイトでみられるはず
The individual AA-ICC corresponds to ICC(A,1) in McGraw and Wong (1996a)
or ICC(2,1) in Shrout and Fleiss (1979). The average AA-ICC corresponds to
ICC(A,k) in McGraw and Wong (1996a) or ICC(2,k) in Shrout and Fleiss (1979).
つまり，ICC(2,1)はindividualの方を読み取ればよいわけです。正確な引用文献も
分かります。
（追記ここまで）*/
こちらのデータにicc23コマンドで，model(3)オプション外して解析してみたところ，
値が一致しました。
これで何となくretestのiccが報告できそうな気がしてきたので，
そのうち，Rの方でもやり方を確かめてみるつもりです。
20150429 individualとaverageの読み取り方追記 たとえば，マルチレベル分析などの際に，職場の部署や
学校のクラス別に，それぞれの合計人数や部署・クラス内ID
を振り分けて，変数として追加したいという場合があります。
ちなみに，同じ状況で層別の平均値や割合については以前
グループごとの平均値や割合を新たな変数として作成
で解説しました。今度は層別の合計人数と層内のIDになります。
この時に便利なのが，_nと_Nです。
_n は，何行目のケースであるか，
_N はケースの総数
をそれぞれ示しています。たとえば，
sysuse auto
で最初の5ケースの車種を表示した場合，
list make in 1/5
+---------------+
| make |
|---------------|
1. | AMC Concord |
2. | AMC Pacer |
3. | AMC Spirit |
4. | Buick Century |
5. | Buick Electra |
+---------------+
となり，この内1番目の AMC Concord だけ表示したい場合，この_nを使うと，
list make if _n==1
+-------------+
| make |
|-------------|
1. | AMC Concord |
+-------------+
となります。一方，_Nの方は，
display _N
74
で，ケースの総数が表示されました。
これらを応用することで，本記事のタイトルにある作業が可能になります。
まず，rep78の各値ごとの合計数をみてみます。 Freq.がその数になります。
tab rep78,mis
Repair |
Record 1978 | Freq. Percent Cum.
------------+-----------------------------------
1 | 2 2.70 2.70
2 | 8 10.81 13.51
3 | 30 40.54 54.05
4 | 18 24.32 78.38
5 | 11 14.86 93.24
. | 5 6.76 100.00
------------+-----------------------------------
Total | 74 100.00
rep78のカテゴリごとにidを振るには，
bysort rep78: gen id_by_rep=_n
次に，rep78のカテゴリごとの人数を変数にするには
bysort rep78: gen n_by_rep=_N
これらをrep78の各カテゴリの最初の5人ずつ表示してみます。
分かりやすいように，sepbyオプションをつけてrep78の値ごとに区切る
ようにします
list rep78 id_by_rep n_by_rep if inlist(id_by_rep,1,2,3,4,5),
sepby(rep78)
+-----------------------------+以前， 層別に平均値と標準偏差（SD）と人数だけを算出 という記事で解説
| rep78 id_by_~p n_by_rep |
|-----------------------------|
1. | 1 1 2 |
2. | 1 2 2 |
|-----------------------------|
3. | 2 1 8 |
4. | 2 2 8 |
5. | 2 3 8 |
6. | 2 4 8 |
7. | 2 5 8 |
|-----------------------------|
11. | 3 1 30 |
12. | 3 2 30 |
13. | 3 3 30 |
14. | 3 4 30 |
15. | 3 5 30 |
|-----------------------------|
41. | 4 1 18 |
42. | 4 2 18 |
43. | 4 3 18 |
44. | 4 4 18 |
45. | 4 5 18 |
|-----------------------------|
59. | 5 1 11 |
60. | 5 2 11 |
61. | 5 3 11 |
62. | 5 4 11 |
63. | 5 5 11 |
|-----------------------------|
70. | . 1 5 |
71. | . 2 5 |
72. | . 3 5 |
73. | . 4 5 |
74. | . 5 5 |
+-----------------------------+
しましたが，層となる変数を組み合わせて，縦一列に平均値だけ欲しいという場合も
よくあります。以下の3通りの方法をまとめます。
sumを使う
tabstatを使う
tableを使う
tabstatを使う
tableを使う
まず，サンプルデータセットを呼び出します。
sysuse auto
■ sumを使う
おなじみのsumを使いbysortとの組み合わせで出す場合，
bysort foreign: sum price mpg weight length
-> foreign = Domestic
Variable | Obs Mean Std. Dev. Min Max
------------+---------------------------------------------
price | 52 6072.423 3097.104 3291 15906
mpg | 52 19.82692 4.743297 12 34
weight | 52 3317.115 695.3637 1800 4840
length | 52 196.1346 20.04605 147 233
---------------------------------------------------------
-> foreign = Foreign
Variable | Obs Mean Std. Dev. Min Max
------------+---------------------------------------------
price | 22 6384.682 2621.915 3748 12990
mpg | 22 24.77273 6.611187 14 41
weight | 22 2315.909 433.0035 1760 3420
length | 22 168.5455 13.68255 142 193
層別変数を2つ組み合わせることもできます。
bysort foreign rep78: sum price mpg weight length
-> foreign = Domestic, rep78 = 1
Variable | Obs Mean Std. Dev. Min Max
------------+---------------------------------------------
price | 2 4564.5 522.5519 4195 4934
mpg | 2 21 4.242641 18 24
weight | 2 3100 523.259 2730 3470
length | 2 189 12.72792 180 198
----------------------------------------------------------
-> foreign = Domestic, rep78 = 2
Variable | Obs Mean Std. Dev. Min Max
------------+---------------------------------------------
price | 8 5967.625 3579.357 3667 14500
mpg | 8 19.125 3.758324 14 24
weight | 8 3353.75 445.9961 2690 3900
length | 8 199.375 13.97894 179 220
(長いので省略，Domesticが終わったら，foreign = Foreignについても出力される)
■ tabstatを使う （20150412 追記）
tabstatは，デフォルトだと下に示すように変数が横に並んで出てきます。
tabstat price mpg weight length, stat(mean) by(foreign)
Summary statistics: mean
by categories of: foreign (Car type)
foreign | price mpg weight length
---------+----------------------------------------
Domestic | 6072.423 19.82692 3317.115 196.1346
Foreign | 6384.682 24.77273 2315.909 168.5455
---------+----------------------------------------
Total | 6165.257 21.2973 3019.459 187.9324
--------------------------------------------------
totalがいらない場合は，nototalオプションをつけます
tabstat price mpg weight length, stat(mean) by(foreign)
nototal
Summary statistics: mean
by categories of: foreign (Car type)
foreign | price mpg weight length
---------+----------------------------------------
Domestic | 6072.423 19.82692 3317.115 196.1346
Foreign | 6384.682 24.77273 2315.909 168.5455
--------------------------------------------------
変数を縦に表示したい場合は，col(stat)オプションをつけます
tabstat price mpg weight length, stat(mean) by(foreign) nototal
col(stat)
foreign | mean
---------+----------
Domestic | 6072.423
| 19.82692
| 3317.115
| 196.1346
---------+----------
Foreign | 6384.682
| 24.77273
| 2315.909
| 168.5455
--------------------
真ん中の線もいらないので，nosepオプションをつけます
tabstat price mpg weight length, stat(mean) by(foreign) nototal col(stat)
nosep
foreign | mean
---------+----------
Domestic | 6072.423
| 19.82692
| 3317.115
| 196.1346
Foreign | 6384.682
| 24.77273
| 2315.909
| 168.5455
--------------------
平均値だけじゃなく，Nも表示したい場合はstat()に追加します
tabstat price mpg weight length, stat(mean
n
) by(foreign) nototal col(stat) nosep
foreign | mean N
---------+--------------------
Domestic | 6072.423 52
| 19.82692 52
| 3317.115 52
| 196.1346 52
Foreign | 6384.682 22
| 24.77273 22
| 2315.909 22
| 168.5455 22
------------------------------
出力には対象となった変数の説明が表上部に
表示されているのですが，表の中に変数名も入れ込みたい場合
はlongstubオプションをつけます
tabstat price mpg weight length, stat(mean) by(foreign) nototal col(stat) nosep
longstub
foreign variable | mean
----------------------+----------
Domestic price | 6072.423
mpg | 19.82692
weight | 3317.115
length | 196.1346
Foreign price | 6384.682
mpg | 24.77273
weight | 2315.909
length | 168.5455
---------------------------------
bysortとの組み合わせもできるので，2要因以上も大丈夫です
bysort foreign:
tabstat price mpg weight length, stat(mean) by(rep78) nototal col(stat) nosep longstub
-> foreign = Domestic
rep78 variable | mean
----------------------+----------
1 price | 4564.5
mpg | 21
weight | 3100
length | 189
2 price | 5967.625
mpg | 19.125
weight | 3353.75
length | 199.375
3 price | 6607.074
mpg | 19
weight | 3442.222
length | 197.8889
4 price | 5881.556
mpg | 18.44444
weight | 3532.222
length | 204.4444
5 price | 4204.5
mpg | 32
weight | 1960
length | 160
---------------------------------
-------------------------------------------------
-> foreign = Foreign
rep78 variable | mean
----------------------+----------
3 price | 4828.667
mpg | 23.33333
weight | 2010
length | 159
4 price | 6261.444
mpg | 24.88889
weight | 2207.778
length | 165.2222
5 price | 6292.667
mpg | 26.33333
weight | 2403.333
length | 172.4444
---------------------------------
■ tableを使う
とにかく平均値だけ欲しいときは，tableを使うと無駄のない出力が出てきます。
縦にそれぞれの変数の平均値だけが並んでいます。
table rep78 foreign , c(mean price mean mpg mean weight mean length )
------------------------------たとえば，100個の変数について性別とのクロス表をいっぺんに出したいという
Repair |
Record | Car type
1978 | Domestic Foreign
----------+-------------------
1 | 4,564.5
| 21
| 3,100
| 189
|
2 | 5,967.6
| 19.125
| 3,353.8
| 199.375
|
3 | 6,607.1 4,828.7
| 19 23.3333
| 3,442.2 2,010
| 197.889 159
|
4 | 5,881.6 6,261.4
| 18.4444 24.8889
| 3,532.2 2,207.8
| 204.444 165.222
|
5 | 4,204.5 6,292.7
| 32 26.3333
| 1,960 2,403.3
| 160 172.444
------------------------------
状況があるとします。以下，Stata12.0で説明します。
性別の変数sexと変数a1 a2 a3 ･･･ a100とあったとき，
クロス表を書くコマンドは，
tab sex a1
tab sex a2
tab sex a3
・
・
・
tab sex a100
tab sex a2
tab sex a3
・
・
・
tab sex a100
と100行にわたってしまいます。でも100行も書いてられないですよね。
実はこれは1行で書くことができます。
for varlist a1-a100: tab sex X
これは，for varlistの後におかれた一連の変数について，その後のXに代入して処理を
繰り返してくれるコマンドです。
a1-a100というのは，連番を意味していて，-(ハイフン)でつながれた間の変数
もすべて含めてくれるので，100個も変数を並べなくて大丈夫です。
もし並べたい場合は，a1 a2 a3 ･･･ a100とスペースを入れて並べます。
このようにfor～ というコマンドは便利なのですが，最近はout of dateだ
そうで，ヘルプにも詳しく書かれておりません。
もしかしたらいずれ使えなくなってしまうかもしれません。
それでは，今後を見据えて最新のコマンドではどう書くようになっているか
というと，最低でも3行必要なように変わっているようです。
foreach x of varlist a1-a100 {
tab sex `x'
}
tab sex `x'
}
または，
foreach x in a1-a100 {
tab sex `x'
}
tab sex `x'
}
です。
これは，foreachの後に置かれた文字や単語（なんでもよい）に一連の変数を
代入してくれるという意味です。for～との違いは，代入先が柔軟になったと
いうところでしょうか。
気をつけなければいけないのは， { と } の位置は必ず上記の通りに
決まっていること，また ` ' という初めての人には見慣れない記号がある
ことです。これらは，
a left quote ( ` )
a right quote ( ' )
という名前のようで，自分の環境だと，半角の状態でshift押しながら@と7キーを押す
と出てきました。
ここでやっていることは，英語では looping というみたいです。
検索すると大体この言葉と一緒にひっかかってきます。
これらの動きが分かってくると，一気にStataへの親しみがわいてきます。
（2015/04/06 追記）
■ 数値のみを繰り返したい時
これは，forvaluesを使います。下記の記事参照。
線形回帰のメモ （いわゆる単回帰分析，重回帰分析）
の，(5)グループごとの回帰（層別）で実行例あり。
（2015/05/07 追記）
また応用編として↓も
連続して変わる数値の繰り返し処理を入れ子で実行する（double loop）
変更履歴
2013/11/17 度数分布表→クロス表に修正
2015/04/06 forvaluesの例追記
2015/05/07 リンク追記，foreachのコマンド例に，クロス表の例として一貫するようsexを追記 縦軸をアウトカム，横軸を時間変数としてアウトカムの推移を1人ずつ
見ていけるグラフを描きます。まずは全体で出してから，治療群別に出したものも
示します。
データセットはBtheBを使います。longデータセット形式の方が
楽なので，まずはデータセットの変換が必要です。用意の仕方は
データセットをwideからlongへ変換する
を参照ください。なお，ここではbdi_preのデータもbdi変数
の中に含めます。なので，変数名の変更の時に，
rename (
bdi_pre
bdi_2m bdi_3m bdi_5m bdi_8m )(bdi0 bdi2 bdi3 bdi5 bdi8 )
とします。
全部で100名いるので，全員グラフにするとかなり時間がかかってしま
います。したがって，まずはidが20以下の20名分（if id<=20）をグラフ化してみます。
縦軸のアウトカムが抑うつ症状bdiで，横軸が時点のmonthです。
xtline bdi if id<=20, t(month) i(id)
次に，これらを1枚のグラフにまとめたものは，
オプションにoverlayをつけると出てきます
xtline bdi if id<=20, t(month) i(id) overlay
さらにオプションに legend(off) を加えると，idと線の色の対応を示して
いる凡例部分がなくなり，グラフが大きくなります。コマンドだけ
書いておきます。
xtline bdi if id<=20, t(month) i(id) legend(off) overlay
治療群別に見たい場合は，ifの後にさらに条件を加えて群ごとに描くことで
出てきます。どちらの群か分かりやすくするため，タイトルも付けます。
xtline bdi if treatment==1 & id<=20, t(month) i(id) title("TAU") legend(off) overlay
xtline bdi if treatment==2 & id<=20, t(month) i(id) title("BtheB") legend(off) overlay
別々に出てきてしまったので，グラフを1枚にまとめたいときは
以下のようにしていきます。まず，グラフをたくさん描いてきたのでいったん消去するために，
graph drop _all
そして，グラフをまとめるためにそれぞれ名前（taとbt）をつけて，合併するコマンド
は以下です。
xtline bdi if treatment==1 & id<=20, t(month) i(id) title("TAU") legend(off) overlay name(ta)
xtline bdi if treatment==2 & id<=20, t(month) i(id) title("BtheB") legend(off) overlay name(bt)
graph combine ta bt
(20150320 修正)
TAUの方のy軸の範囲がBtheBとそろってないのを修正
xtline bdi if treatment==1 & id<=20, t(month) i(id) title("TAU") legend(off) overlay name(ta)
ylabel(0(10)50)
xtline bdi if treatment==2 & id<=20, t(month) i(id) title("BtheB") legend(off) overlay name(bt)
graph combine ta bt
結論から言うとtabでの解決策は不明です。ただ，いくつか
代替策があったのでメモしておきます。
自分でデータを入力して試してみます。まず
doedit
でdoファイルエディタを開き，以下を書いて実行します。
input str60 moji
"あいう（全角１０字）"
"あいうえおかきく（全角１５字）"
"あいうえおかきくけこさしす（全角２０字）"
"あいうえおかきくけこさしすせそたちつ（全角２５字）"
"あいう（全角１０字）"
"あいう（全角１０字）"
"あいう（全角１０字）"
"あいうえおかきく（全角１５字）"
"あいうえおかきく（全角１５字）"
"あいうえおかきくけこさしす（全角２０字）"
"."
"."
end
"あいう（全角１０字）"
"あいうえおかきく（全角１５字）"
"あいうえおかきくけこさしす（全角２０字）"
"あいうえおかきくけこさしすせそたちつ（全角２５字）"
"あいう（全角１０字）"
"あいう（全角１０字）"
"あいう（全角１０字）"
"あいうえおかきく（全角１５字）"
"あいうえおかきく（全角１５字）"
"あいうえおかきくけこさしす（全角２０字）"
"."
"."
end
str60というのは，文字系変数で半角60字までの長さだという指定
です。詳しくは
help data types
それでは，データセットを確認してみます。
list
+--------------------------------------------------+
| moji |
|--------------------------------------------------|
1. | あいう（全角１０字） |
2. | あいうえおかきく（全角１５字） |
3. | あいうえおかきくけこさしす（全角２０字） |
4. | あいうえおかきくけこさしすせたちつ（全角２０字） |
5. | あいう（全角１０字） |
|--------------------------------------------------|
6. | あいう（全角１０字） |
7. | あいう（全角１０字） |
8. | あいうえおかきく（全角１５字） |
9. | あいうえおかきく（全角１５字） |
10. | あいうえおかきくけこさしす（全角２０字） |
|--------------------------------------------------|
11. | . |
12. | . |
+--------------------------------------------------+
ここでは文字が長くても表示されています。
tab moji,mis
moji | Freq. Percent Cum.
----------------------------------------+-----------------------------------
. | 2 16.67 16.67
あいう（全角１０字） | 4 33.33 50.00
あいうえおかきく（全角１５字） | 3 25.00 75.00
あいうえおかきくけこさしす（全角２ ０ ・. | 2 16.67 91.67
あいうえおかきくけこさしすせそたち つ ・. | 1 8.33 100.00
----------------------------------------+-----------------------------------
Total | 12 100.00
この赤字の部分のように，度数を確認した時に，長い文字列がとぎれてしまうのを
何とかしたいわけです。全角18字までが限界なようです。
もう少しだけ長い文字列が表示できるのが，tableです。
table moji, mis
-----------------------------------------------------
moji | Freq.
-----------------------------------------+-----------
. | 2
あいう（全角１０字） | 4
あいうえおかきく（全角１５字） | 3
あいうえおかきくけこさしす（全角２０字） | 2
あいうえおかきくけこさしすせそたちつ（ 全 | 1
-----------------------------------------------------
こちらは赤字部分までの全角20字（半角40字）までは表示されましたが
それが限界のようです。
このように，度数は分かっているので，この文字列が判読できてコピーできれば
よいわけです。そこで，苦肉の策として，listを工夫してみます。
doファイルエディタの方で，重複する値を一時的にdropして，
listを出したら元に戻す方法でやってみます。それぞれの値が1ケースだけ
残るというわけです。
preserve
duplicates drop moji, force
list, sep(0)
restore
duplicates drop moji, force
list, sep(0)
restore
+----------------------------------------------------+
| moji |
|----------------------------------------------------|
1. | あいう（全角１０字） |
2. | あいうえおかきく（全角１５字） |
3. | あいうえおかきくけこさしす（全角２０字） |
4. | あいうえおかきくけこさしすせそたちつ（全角２５字） |
5. | . |
+----------------------------------------------------+
欠損の位置が逆になりましたが，後は同じ順番で並んでいます。listを出す際に
重複を削りましたが，preserve～restoreではさんで一時的に実行して元に戻して
いるので，データセットの中身は変わっていません。tabやtableと並び順が変わっても，
Excelとかにはり付けてソートすればそろいます。 まず以下のコマンドを打ち込み，例となるデータセットを
作成します。データセットは，1日の内の出勤時刻（h1とm1）と
退勤時刻（h2とm2）の計4つの変数があるイメージです。
時と分が異なる列に入力されているケースです。
使用できる値の確認のため，24時を超えた26時という値も入れています。
doファイルエディタを開いて，以下のコマンドを入力し実行します。
input h1 m1 h2 m2
9 0 17 30
9 0 17 0
9 0 26 0
end
list
+-------------------------+
| h1 m1 h2 m2 |
|---------------------------|
1. | 9 0 17 30 |
2. | 9 0 17 0 |
3. | 9 0 26 0 |
+-------------------------+
それぞれ時刻の変数を作って差（つまり勤務時間）を求める
方法をいろいろ試してみたらうまくいかなかったので，単純に
分単位に揃えて分の差を算出する方法を紹介します。
新たにin_minとout_minの2つの変数を作成します。
h1，h2は時刻の変数なので，60を掛けることで分を表します。
m1とm2は元々分なのでそのまま足すだけです。
gen in_min=60*h1+m1
gen out_min=60*h2+m2
gen out_min=60*h2+m2
それぞれ時刻を分換算したので，あとは差を求めるだけです。
difminという変数名にしています。
gen difmin= out_min - in_min
list
list
+-------------------------------------------------------------------+
| h1 m1 h2 m2 in_min out_min difmin |
|-------------------------------------------------------------------|
1. | 9 0 17 30 540 1050 510 |
2. | 9 0 17 0 540 1020 480 |
3. | 9 0 26 0 540 1560 1020 |
+------------------------------------------------------------------+
1番目の方は8時間半の差なので，分換算したら，60*8+30=510
なので，これで時間差difminが分で求まりました。
ちなみに，3番目の方は，17時間の差なので，分換算
すると60*17=1020です。このように時刻が24時を超えても
経過時間が反映する変数となっていれば問題ありません。
時間差の平均値も簡単に求まります。
sum difmin
Variable | Obs Mean Std. Dev. Min Max
----------------+-------------------------------------------------------------------------
difmin | 3 670 303.4798 480 1020
平均670分なので，11時間10分が平均と解釈できます。
[変更履歴]
2015/02/06 データ入力部分と listの後に変数名1つだけだった部分修正 ヒストグラムの描き方での基本メモです。
autoデータを使用します。
sysuse auto
(1)カテゴリデータの場合
変数はrep78です。
hist rep78
デフォルトのヒストグラムだと，このように見やすくない感じです。
少しずつ改良していくには，オプションを足していきます。
縦軸が密度になっているので，まずは縦軸を度数にしてみます。
hist rep78, frequency
frequencyはfreqと略しても同様ですので，以下は
freqを使います。
次に，discreteオプションを使いデータが離散変数であると
いう指定をします。要は，横軸の1単位ごとにグラフを
表示したい場合に使います。
hist rep78, freq discrete
discreteも同様にdiscと略せるので，以下略して使います。
これでかなり見やすくなりました。最後に，
それぞれの人数をグラフ上に示してます。
hist rep78, freq disc addlabel
(2)連続量の場合
変数はmpgです。
hist mpg, freq
横軸の目盛りをもっと詳細にしたい場合は，まず
最小値と最大値を確認するために
sum mpg
Variable | Obs Mean Std. Dev. Min Max
----------------+-----------------------------------------------------------------------
mpg | 74 21.2973 5.785503 12 41
この情報を元に，横軸を5単位ずつ表示させるためには，xlabelオプション
を使います。xlabelの（ ）内に，横軸の端から端までの数値を入れその間に
さらに（ ）を入れて，何単位ごとの表示か指定するという構造です。
せっかく目盛りが細かくなるので，ヒストグラムも横軸の1単位ずつ表示させる
ためdiscオプションもつけます。
hist mpg, freq xlabel(12(5)41) disc
次は，foreignの国産／外国車でヒストグラムを層別に
描いてみます
hist mpg, freq by(foreign)
同じグラフを今度は上下に並べるように描きます。
hist mpg, by(foreign, col(1))
(3)ヒストグラムの重ね描き
いつも参照しているUCLAの統計コンサルのサイトより
How can I overlay two histograms?
ただしここではデータをautoデータにしてやってみます。
まず，foreignのラベル内訳と数字を確認するため，
tab foreign
Car type | Freq. Percent Cum.
----------------+-----------------------------------
Domestic | 52 70.27 70.27
Foreign | 22 29.73 100.00
----------------+-----------------------------------
Total | 74 100.00
tab foreign, nolabel
Car type | Freq. Percent Cum.
---------------+-----------------------------------
0 | 52 70.27 70.27
1 | 22 29.73 100.00
---------------+-----------------------------------
Total | 74 100.00
したがって，0:Domestic, 1:Foreign という対応です。
デフォルトのまま描くと同じ色になってしまうので，Domesticの方を緑にします。
また，数行に渡るので，doファイルエディタを使います。legendオプションは
凡例です。これをつけないと区別がつきにくくなるのでつけています。
twoway (hist mpg if foreign==0, freq disc color(green))
///
(hist mpg if foreign==1, freq disc ///
legend(order(1 " Domestic" 2 "Foreign" )))
(hist mpg if foreign==1, freq disc ///
legend(order(1 " Domestic" 2 "Foreign" )))
これだと重なってしまって分からない部分があるので，
fcolor(none)とlcolor(black)のオプションを使って描きます。
twoway (hist mpg if foreign==0, freq disc color(green))
///
(hist mpg if foreign==1, freq disc fcolor(none) lcolor(black) ///
legend(order(1 " Domestic" 2 "Foreign" )))
(hist mpg if foreign==1, freq disc fcolor(none) lcolor(black) ///
legend(order(1 " Domestic" 2 "Foreign" )))
階層的重回帰とは，一連の説明変数を段階を踏んで投入していき
R二乗の変化量を出して検定を行う，といった感じの分析を指します。
自分はあまり使わないのですが，聞かれたのでやり方を調べて
みました。
autoデータセットを使います。Stataのバージョンは13.1です。
sysuse auto
単純に階層的に説明変数を追加していく回帰モデルを
並べると，
reg mpg price
reg mpg price foreign weight
reg mpg price foreign weight length
reg mpg price foreign weight
reg mpg price foreign weight length
のようになります。こうした場合，R二乗の変化量の検定を簡単に
やってくれるコマンドがあります。
nestreg: reg mpg price (foreign weight) (length), nohead
最終的に用いる説明変数全部を入れたregコマンドについて，
投入する一連の変数（ブロック）ごとに( )でくくり，regの前に
nestreg: というのをつけるだけです。出力で分散分析表を出すと長く
なるので，noheadオプションで抑制しています。
Block 1: price
---------------------------------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+-----------------------------------------------------------------------------------------
price | -.0009192 .0002042 -4.50 0.000 -.0013263 -.0005121
_cons | 26.96417 1.393952 19.34 0.000 24.18538 29.74297
---------------------------------------------------------------------------------------------------------
Block 2: foreign weight
----------------------------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+--------------------------------------------------------------------------------------
price | .0000566 .0001922 0.29 0.769 -.0003268 .00044
foreign | -1.855891 1.289063 -1.44 0.154 -4.426846 .7150641
weight | -.0067758 .0009048 -7.49 0.000 -.0085805 -.0049712
_cons | 41.95948 2.377726 17.65 0.000 37.21725 46.7017
------------------------------------------------------------------------------------------------------
Block 3: length
--------------------------------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------------------------------
price | -.0000374 .0002009 -0.19 0.853 -.0004381 .0003634
foreign | -1.574443 1.292234 -1.22 0.227 -4.15238 1.003494
weight | -.0041499 .0019865 -2.09 0.040 -.0081128 -.0001871
length | -.086156 .058149 -1.48 0.143 -.2021601 .029848
_cons | 50.71772 6.364007 7.97 0.000 38.02187 63.41357
------------------------------------------------------------------------------------------------------
+--------------------------------------------------------------------------------+
| | Block Residual Change |
| Block | F df df Pr > F R2 in R2 |
|----------+----------------------------------------------------------------------|
| 1 | 20.26 1 72 0.0000 0.2196 |
| 2 | 46.08 2 70 0.0000 0.6631 0.4435 |
| 3 | 2.20 1 69 0.1430 0.6735 0.0104 |
+--------------------------------------------------------------------------------+
カテゴリ変数は非対応みたいですが，2値変数ならOKなのでforeignが入っています。
標準編回帰係数で出したい場合は，さらにオプションにbetaを入れます。
検証的（確証的／確認的）因子分析：confirmatory factor analysis(CFA)
の続きです。
因子間相関の有無などをどう指定すればよいかについて説明します。
以下の話を理解するには，SEMの用語の
内生変数（endogenous variable）と外生変数（exogenous variable）
という言葉を知っておく必要があります。簡単に言うと，内生変数は，
矢印を受ける変数で，外生変数は矢印を受けておらず，発する側である
変数を指します。
StataのsemコマンドによるCFAでは，何も指定しない場合，デフォルトで
因子間相関が仮定されています。これを，無相関にしたい場合は，
共分散構造を指定する必要があり，次のオプションを追加します。
covstruct(_LEx, diagonal)
前回のコマンドに加えるだけなので，因子間相関を仮定しない場合は
下記のようになります。
sem (A1-A5 <- Agree)(C1-C5 <- Conscientious)(E1-E5 <- Extraversion)
///
(N1-N5 <- Neuroticism)(O1-O5 <- Openness), nocapslatent ///
latent(Agree Conscientious Extraversion Neuroticism Openness) covstruct(_LEx, diagonal)
(N1-N5 <- Neuroticism)(O1-O5 <- Openness), nocapslatent ///
latent(Agree Conscientious Extraversion Neuroticism Openness) covstruct(_LEx, diagonal)
_LEx とは，全ての潜在外生変数を指すので，ここでは因子のことです。
_LEx の代わりに，ここに因子名をすべて書いて
covstruct(Agree Conscientious Extraversion Neuroticism Openness, diagonal)
としてもよいのですが，長くなるので，_LExで省略しているというわけです。
diagonalは，全ての分散は制約がなく，全ての共分散は0に固定という指定です。
この辺りの事は詳しくは
help sem and gsem option covstructure
とコマンドを打ってヘルプを確認すると分かります。
デフォルトでは，
semの外生変数の共分散構造は，covstruct(_Ex, unstructured)，
semとgsem共に，誤差変数の共分散構造は，covstruct(e._En, diagonal)
となっているそうです。 Stata 13.1での検証的（確証的／確認的）因子分析について，共分散構造分析
を用いた分析の仕方をまとめます。主には
help two-factor measurement model
のコマンドで参照されるマニュアルの記述を参考にしています。
データは，現実に扱うデータのタイプによく合致していた
Rのパッケージpsychのbfiデータを使います。
RのデータセットをStataデータセットに変換するには↓参照
EZRでパッケージ中のRデータセットをStataデータセットに変換する
■bfiデータの説明
まず，Rにpsychパッケージが入っていることが前提ですが，
Rの方で
library(psych)
help(bfi)
と書いて実行すると説明が読めます。
25 Personality items representing 5 factors
5因子を表すパーソナリティに関する25項目
とあり，
国際パーソナリティ項目プールからの2800名分のデータ。
それぞれの因子に対応した質問項目A1からA5，C1からC5，E1からE5，N1からN5，O1からO5
の変数に加え，性，教育歴，年齢の変数もあり
といった内容になっています。
想定される因子は，5項目ずつにあり
Agree →A1からA5
Conscientious →C1からC5
Extraversion →E1からE5
Neuroticism →N1からN5
Openness →O1からO5
Conscientious →C1からC5
Extraversion →E1からE5
Neuroticism →N1からN5
Openness →O1からO5
となっていて，因子名の頭文字が変数名と対応しています。
■解析手順
ここからはStataでのコマンドです。
まず，Stataデータセットに変換したbfiデータを開きます。
次に，記述統計を確認します。
sum A1-O5
とハイフンを使えば簡単に全項目の記述統計が見られます。
途中まで示してみると，
Variable | Obs Mean Std. Dev. Min Max
--------------+------------------------------------------------------------------------
A1 | 2784 2.413434 1.407737 1 6
A2 | 2773 4.80238 1.17202 1 6
A3 | 2774 4.603821 1.301834 1 6
A4 | 2781 4.699748 1.479633 1 6
A5 | 2784 4.560345 1.258512 1 6
-------------+--------------------------------------------------------
C1 | 2779 4.502339 1.241347 1 6
C2 | 2776 4.369957 1.318347 1 6
C3 | 2780 4.303957 1.288552 1 6
C4 | 2774 2.553353 1.375118 1 6
C5 | 2784 3.296695 1.628542 1 6
-------------+-------------------------------------------------------------------------
という感じです。
Obsがばらばらなので，適度に欠損値がみられるデータということになります。
欠損のパターンをみたい場合は，↓参照
データの欠損を数える
さて，いよいよCFAの分析に入ります。
SEM Builderを使ってお絵描き風にやってもよいのですが，
コマンドで書いた方がはるかに楽なので，次のようにします。
途中で改行してるので，do-file エディタに以下のように書きます。
エディタはメニューのアイコン押すか，doeditというコマンドで開きます
sem (A1-A5 <- Agree)(C1-C5 <- Conscientious)(E1-E5 <- Extraversion)
///
(N1-N5 <- Neuroticism)(O1-O5 <- Openness), nocapslatent ///
latent(Agree Conscientious Extraversion Neuroticism Openness)
(N1-N5 <- Neuroticism)(O1-O5 <- Openness), nocapslatent ///
latent(Agree Conscientious Extraversion Neuroticism Openness)
簡単に解説すると，(A1-A5 <-Agree)というのは，A1-A5の項目にAgreeという
潜在変数（因子）が対応しているという意味で，それが5因子分繰り返されています。
通常は，潜在変数の頭文字が大文字，観測変数の頭文字が小文字であることを
想定しているので，変数の頭文字が大文字だったりする場合，nocapslatent
オプションとlatentオプションをつけるみたいです。latent（ ）の中に因子名が
並んでいますが，ここで潜在変数を明示しています。
結果は長いので省きますが，係数の一覧等が出てくるはずです。
標準化係数を出すには，次のコマンドを実行します。
sem, standardized
次に，適合度指標一覧を算出します。
estat gof, stats(all)
----------------------------------------------------------------------------
Fit statistic | Value Description
-------------------------------+------------------------------------------------------
Likelihood ratio |
chi2_ms(265) | 4165.467 model vs. saturated
p > chi2 | 0.000
chi2_bs(300) | 18222.116 baseline vs. saturated
p > chi2 | 0.000
-------------------------------+------------------------------------------------------
Population error |
RMSEA | 0.078 Root mean squared error of approximation
90% CI, lower bound | 0.076
upper bound | 0.080
pclose | 0.000 Probability RMSEA <= 0.05
-------------------------------+------------------------------------------------------
Information criteria |
AIC | 199850.476 Akaike's information criterion
BIC | 200343.316 Bayesian information criterion
-------------------------------+------------------------------------------------------
Baseline comparison |
CFI | 0.782 Comparative fit index
TLI | 0.754 Tucker-Lewis index
-------------------------------+------------------------------------------------------
Size of residuals |
SRMR | 0.073 Standardized root mean squared residual
CD | 0.999 Coefficient of determination
-----------------------------------------------------------------------------------------
一般的には，カイ二乗値，CFI, TLI, RMSEAとその90%信頼区間，SRMRなどが分かれば
よいので，これで論文に記載する情報がとりあえずそろいます。
因子間相関については
検証的（確証的／確認的）因子分析 (2)因子間相関
全体的にUCLAの統計コンサルページの
How can I understand a categorical by continuous interaction? (Stata 12)
を参考にしています。
まずはautoデータセットを開きます。
sysuse auto
車の値段に対する，重量と生産国（国産 vs 外国産）の影響と
その交互作用を解析します。
reg price c.weight##foreign, noheader
ここで，c.weightとしてc.をつけているのは，weightを
連続量だと定義するためです。noheaderは，分散分析表などを
出力しないオプションです。以下のような出力が出てきます。
-----------------------------------------------------------------------------------------------------
price | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-----------------+-----------------------------------------------------------------------------------
weight | 2.994814 .4163132 7.19 0.000 2.164503 3.825124
|
foreign |
Foreign | -2171.597 2829.409 -0.77 0.445 -7814.676 3471.482
|
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
-----------------------------------------------------------------------------------------------------
交互作用項foreign#c.weightが有意でしたので，その関係性を詳しく
確認したいとします。weightの一定の区分ごとに，予測される
値段をグラフにしてみます。まず，weightの記述統計を確認すると，
sum weight
Variable | Obs Mean Std. Dev. Min Max
----------------+------------------------------------------------------------------------
weight | 74 3019.459 777.1936 1760 4840
と値の範囲が分かりますので，とりあえず重さ500ごとに
区分してみましょう。そして，国産 vs 外国産で重量の
段階ごとに予測される値段（Margin:調整済み平均）は以下の
コマンドで出力されます。
margins foreign, at(weight=(1760(500)4840))
//weightの後の数値の意味は，最小値（500ごとに区切り）最大値，という意味です。
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を用いてグラフ化します。
marginsplot
グラフを見ると，重量が大きくなるに従い，外国産の方がより値段の上がり方
が大きくなる関係が見られました。ただ，信頼区間の幅も大きくなっているので
ばらつきが大きいことには注意が必要です。
これに散布図も組み合わせると，よりデータの把握に役立ちます。それには
以下のようにオプションをつけます。
marginsplot, addplot(scatter price weight)
marginsの解釈については，
How can I graph the results of the margins command? (Stata 12)
にもあります。この解説によると，
SASではleast squares means (lsmeans)とよばれ，
SPSSでは estimated marginal means (emmeans) とよばれる
と説明があります。
介入研究の解析での線形混合モデル（1） の続きです。
解析後のデータを視覚化して治療効果の違いを把握します。
交互作用が有意でないこともあり，明確な違いを示すグラフには
なりませんが，元の論文でも，推定結果を元にしたらしいmean profilesを
示したグラフを示していたので，こちらでも似たようなものを描いて
みます（preの値を入れてない点は明らかに違いますが）。
mixedを実行した後に，群×時間のそれぞれの周辺効果（調整平均）が
以下のコマンドを打つことで算出できます。
margins treatment#month
----------------------------------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
--------------------+------------------------------------------------------------------------------------
treatment#month |
TAU#2 | 19.46667 1.619558 12.02 0.000 16.29239 22.64094
TAU#3 | 17.85983 1.693887 10.54 0.000 14.53987 21.17978
TAU#5 | 16.25255 1.763116 9.22 0.000 12.7969 19.70819
TAU#8 | 13.52009 1.81302 7.46 0.000 9.966633 17.07354
BtheB#2 | 14.71154 1.506611 9.76 0.000 11.75864 17.66444
BtheB#3 | 13.69984 1.617249 8.47 0.000 10.53009 16.86959
BtheB#5 | 12.83243 1.697532 7.56 0.000 9.505326 16.15953
BtheB#8 | 12.09091 1.721696 7.02 0.000 8.716448 15.46537
---------------------------------------------------------------------------------
これをグラフ化すると以下のように治療効果の違いが視覚的に
分かるようになります。
marginsplot, x(month)
なんだか，TAUの方が下がり方が急な印象。8か月後は群間にあまり
違いがなさそうです。 autoデータセットを使って説明します。
sysuse auto
変数rep78の内，1と3と5を1とした変数rep78rを作りたいとします
(１) genとreplaceを使う場合
gen rep78r=0
replace rep78r=1 if (rep78==1 | rep78==3 | rep78==5)
replace rep78r=1 if (rep78==1 | rep78==3 | rep78==5)
と書けばよいです *1 。結果は以下の通り
tab rep78 rep78r
Repair |
Record | rep78r
1978 | 0 1 | Total
------------+-----------------------+----------
1 | 0 2 | 2
2 | 8 0 | 8
3 | 0 30 | 30
4 | 18 0 | 18
5 | 0 11 | 11
------------+------------------------+----------
Total | 26 43 | 69
この場合は1と3と5の3つしかないので大変ではないですが，
選びたい値がたくさんあるとき，いちいち
rep78==1 | rep78==3・・・と繰り返すのは大変です。
できれば，重複している所は書かずに，選びたい値だけを
指定したいと思います。そんな時に使えるのがinlist( )関数です。
inlist(変数名, 値1, 値2...) と書いていきます。
gen rep78r2=0
replace rep78r2=1 if inlist(rep78,1,3,5)
replace rep78r2=1 if inlist(rep78,1,3,5)
結果はもちろん上と同様のものになります。
tab rep78 rep78r2
Repair |
Record | rep78r2
1978 | 0 1 | Total
--------------+-----------------------+----------
1 | 0 2 | 2
2 | 8 0 | 8
3 | 0 30 | 30
4 | 18 0 | 18
5 | 0 11 | 11
--------------+-----------------------+----------
Total | 26 43 | 69
ただし，inlist( ) の中に249個以上入れるとエラーに
なるみたいなので，長い場合は分割してreplaceを複数繰り返す
必要があります。
参照: Functions: inlist()
上で説明したやり方は，データが文字型の変数でも使えます。
makeという変数は，メーカーと車種の情報が入っているので，
日本車っぽい名前が入っているのを手動で抜きだしてみました
"Honda Accord"
"Honda Civic"
"Mazda GLC"
"Subaru"
"Toyota Celica"
"Toyota Corolla"
"Toyota Corona"
これらのみ1,他は0とした変数japanを作りたい場合は以下のようにします。
長いので，途中で改行してDo-file editorの方で実行します。
gen japan=0
replace japan=1 if inlist(make,"Honda Accord","Honda Civic", ///
"Mazda GLC","Subaru","Toyota Celica","Toyota Corolla","Toyota Corona")
replace japan=1 if inlist(make,"Honda Accord","Honda Civic", ///
"Mazda GLC","Subaru","Toyota Celica","Toyota Corolla","Toyota Corona")
(2) recodeを使う場合
recode rep78 (1 3 5=1) (2 4=0), gen(rep78r)
または
recode rep78 (1 3 5=1) (else=0) if rep78~=., gen(rep78r)
でも同じことができます。recodeだと( )内はたくさん指定しても
大丈夫みたいでした。
注：
*1 このままだと欠損もrep78rで0になってしまうので，これを
防ぎたい場合はgen rep78r=0 if rep78~=.とする。
【履歴】
2014/07/07 記事作成
2014/07/21 文字型の扱い追記 ユーザー作成コマンドのci2をインストールすると，
相関係数の信頼区間を算出できます。
インターネットにつながっている状態で，
findit ci2
と打ち，sg159_1というリンクをクリックするとインストール
画面になります。
Stata 13.1では問題なく動きました。
使い方を，autoデータセットで確認します。
sysuse auto
まずはmpgとweightのピアソン積率相関係数を普通に出すと
次のようになります。
pwcorr mpg weight
mpg weight
mpg 1.0000
weight -0.8072 1.0000
そして，信頼区間を出すには以下のように打ちます
ci2 mpg weight, corr
Confidence interval for Pearson's product-moment correlation
of mpg and weight, based on Fisher's transformation.
Correlation = -0.807 on 74 observations (95% CI: -0.874 to -0.710)
オプションにspearmanをつければ，スピアマンの順位相関にも対応します。
また，このコマンドに含まれているcii2の方を使うと，相関係数と人数が
分かればその数値だけでも信頼区間が算出できます。
詳しくは
help ci2
年代別など，層別に平均値と標準偏差を出したい場面は多いです。
サンプルデータを使いながら説明します。
sysuse auto
mpgの平均値，標準偏差，人数をforeign別に出したい
とします。Stata（バージョン13.1）では3つのやり方があります。
【算出したい変数が1つのとき】
以下の3つでどれでも同じ結果が得られます。
tabulate foreign, sum(mpg)
table foreign, c(mean mpg sd mpg n mpg)
tabstat mpg, stat(mean sd n) by(foreign)
table foreign, c(mean mpg sd mpg n mpg)
tabstat mpg, stat(mean sd n) by(foreign)
それぞれ特徴などを説明していきます。
最もシンプルなものはtabulate（略してtabでOK）です。
tab foreign, sum(mpg)
| Summary of Mileage (mpg)
Car type | Mean Std. Dev. Freq.
-----------------+-------------------------------------------------
Domestic | 19.826923 4.7432972 52
Foreign | 24.772727 6.6111869 22
-----------------+-------------------------------------------------
Total | 21.297297 5.7855032 74
より柔軟に出すにはtableを使います。オプションのc( )の所
に，「出したい統計量 変数名」のセットをいくつも書いていく形式です。
table foreign, c(mean mpg sd mpg n mpg)
-------------------------------------------------------------------
Car type | mean(mpg) sd(mpg) N(mpg)
---------------+--------------------------------------------------
Domestic | 19.8269 4.743297 52
Foreign | 24.7727 6.611187 22
---------------------------------------------------------------------
同じく柔軟なのはtabstatです
tabstat mpg, stat(mean sd n) by(foreign)
Summary for variables: mpg
by categories of: foreign (Car type)
foreign | mean sd N
--------------+-------------------------------------------
Domestic | 19.82692 4.743297 52
Foreign | 24.77273 6.611187 22
--------------+------------------------------------------
Total | 21.2973 5.785503 74
----------------------------------------
それぞれヘルプをみるともっと細かい設定について解説されています。
【算出したい変数が複数あるとき】
複数の変数について出したいときは，tabstatが便利です。
4つの変数mpg price weight lengthそれぞれの平均，SD，人数を
foreign別に出したいときは
tabstat mpg price weight length, stat( mean sd n ) by(foreign)
Summary statistics: mean, sd, N
by categories of: foreign (Car type)
foreign | mpg price weight length
--------------+-----------------------------------------------------------
Domestic | 19.82692 6072.423 3317.115 196.1346
| 4.743297 3097.104 695.3637 20.04605
| 52 52 52 52
--------------+--------------------------------------------------------------
Foreign | 24.77273 6384.682 2315.909 168.5455
| 6.611187 2621.915 433.0035 13.68255
| 22 22 22 22
--------------+----------------------------------------
Total | 21.2973 6165.257 3019.459 187.9324
| 5.785503 2949.496 777.1936 22.26634
| 74 74 74 74
--------------------------------------------------------------------------------
表の上の所に，算出した統計量が表示されています。
繰り返しの命令を使って，tabulateやtableでもできます。
foreach x of varlist mpg price weight length {
tab foreign, sum(`x')
}
tab foreign, sum(`x')
}
| Summary of Mileage (mpg)
Car type | Mean Std. Dev. Freq.
-----------------+------------------------------------------------
Domestic | 19.826923 4.7432972 52
Foreign | 24.772727 6.6111869 22
-----------------+------------------------------------------------
Total | 21.297297 5.7855032 74
| Summary of Price
Car type | Mean Std. Dev. Freq.
-----------------+------------------------------------
Domestic | 6,072.423 3,097.104 52
Foreign | 6,384.682 2,621.915 22
-----------------+------------------------------------
Total | 6,165.257 2,949.496 74
| Summary of Weight (lbs.)
Car type | Mean Std. Dev. Freq.
-----------------+------------------------------------------------
Domestic | 3,317.115 695.36374 52
Foreign | 2,315.909 433.00345 22
-----------------+------------------------------------------------
Total | 3,019.459 777.19357 74
| Summary of Length (in.)
Car type | Mean Std. Dev. Freq.
-----------------+-------------------------------------------------
Domestic | 196.13462 20.046054 52
Foreign | 168.54545 13.682548 22
-----------------+-------------------------------------------------
Total | 187.93243 22.26634 74
人数やtotalの部分も不要なら，tableを使ってもっと
シンプルな出力ができます
foreach x of varlist mpg price weight length {
table foreign, c(mean `x' sd `x' )
}
table foreign, c(mean `x' sd `x' )
}
--------------------------------------------------
Car type | mean(mpg) sd(mpg)
---------------+-----------------------------------
Domestic | 19.8269 4.743297
Foreign | 24.7727 6.611187
----------------------------------------------------
-------------------------------------------------
Car type | mean(price) sd(price)
---------------+----------------------------------
Domestic | 6,072.4 3097.104
Foreign | 6,384.7 2621.915
---------------------------------------------------
-------------------------------------------------------
Car type | mean(weight) sd(weight)
---------------+--------------------------------------
Domestic | 3,317.1 695.3638
Foreign | 2,315.9 433.0034
-----------------------------------------------------
---------------------------------------------------
Car type | mean(length) sd(length)
---------------+------------------------------------
Domestic | 196.135 20.04605
Foreign | 168.545 13.68255
-----------------------------------------------------
（20150327 追記）
こちらのやり方も参照
層別に複数変数の平均値を縦一列に算出
データの欠損を数える
で解説したように，縦断データは脱落による欠損が多く発生する
ため，すべての時点でデータがある者に限定した古典的な反復測定ANOVA等
では，解析の人数が少なくなってしまう問題がありました。
そこで，よく行われてきた方法は，前の時点と同じ値を次の時点に
代入するlast observation carried forward (LOCF)です。
ただしバイアスが問題となるので，紹介するのはあくまで過去に行われてきた
方法をやってみるためだけのものです。
Stataでは，これを簡単にやってくれるユーザー作成コマンドがあります。
本記事ではStata 13.1を使用してます。インストールは次のコマンドで
findit carryforward
それでは， データの欠損を数える で使用した同じデータセットBtheB
で試してみます。
まず，全時点が欠損である者が3名いたので，これらはこの後の解析
から除外しておきます。
drop if bdi2==. & bdi3==. & bdi5==. & bdi8==.
次に，
データセットをwideからlongへ変換する
でやったように，wideからlong形式に変換します。
準備ができたら，以下のコマンドでLOCFを行います。
ここでは，代入された値を持つ変数，bdinを新しく作成します。
carryforward bdi, gen(bdin)
確認してみます。4行で1名分のデータなので，分割線を4つずつの
指定にします。
list in 1/20, sep(4)
+----------------------------------------------------------------------------+
| id month drug length treatm~t bdi_pre bdi bdin |
|------------------------------------------------------------------------------|
1. | 1 2 No >6m TAU 29 2 2 |
2. | 1 3 No >6m TAU 29 2 2 |
3. | 1 5 No >6m TAU 29 . 2 |
4. | 1 8 No >6m TAU 29 . 2 |
|------------------------------------------------------------------------------|
5. | 2 2 Yes >6m BtheB 32 16 16 |
6. | 2 3 Yes >6m BtheB 32 24 24 |
7. | 2 5 Yes >6m BtheB 32 17 17 |
8. | 2 8 Yes >6m BtheB 32 20 20 |
|------------------------------------------------------------------------------|
9. | 3 2 Yes <6m TAU 25 20 20 |
10. | 3 3 Yes <6m TAU 25 . 20 |
11. | 3 5 Yes <6m TAU 25 . 20 |
12. | 3 8 Yes <6m TAU 25 . 20 |
|------------------------------------------------------------------------------|
13. | 4 2 No >6m BtheB 21 17 17 |
14. | 4 3 No >6m BtheB 21 16 16 |
15. | 4 5 No >6m BtheB 21 10 10 |
16. | 4 8 No >6m BtheB 21 9 9 |
|------------------------------------------------------------------------------|
17. | 5 2 Yes >6m BtheB 26 23 23 |
18. | 5 3 Yes >6m BtheB 26 . 23 |
19. | 5 5 Yes >6m BtheB 26 . 23 |
20. | 5 8 Yes >6m BtheB 26 . 23 |
+-----------------------------------------------------------------------------+
idの1,3,5のbdinの列をみると，bdiで欠損があった時点について，
前の値が代入されていることが確認できます。 ここでは，介入研究のデータセットで解説します。
使用するデータセットとして，
データセットをwideからlongへ変換する
で解説した通り，RのHSAUR2パッケージのデータセットBtheBを
使ってやってみます。Stataのバージョンは13.1です。
データをよみこんだら，準備としてIDづけとlongへの変換用に
変数名の変更は次の通り
egen id=seq()
rename ( bdi_2m bdi_3m bdi_5m bdi_8m )( bdi2 bdi3 bdi5 bdi8 )
rename ( bdi_2m bdi_3m bdi_5m bdi_8m )( bdi2 bdi3 bdi5 bdi8 )
この時点では，縦断データがwide形式，つまりそれぞれ別の変数
として入っている形です。
さて，まずBDIの値の2か月後，3か月後，5か月後，8か月後の
欠損数を一覧にしてみます。次のコマンドで，欠損のある変数の
結果のみ表示されます。
misstable sum
Obs<.
+-------------------------
| | Unique
Variable | Obs=. Obs>. Obs<. | values Min Max
------------+-------------------------------------+-------------------------
bdi2 | 3 97 | 37 0 48
bdi3 | 27 73 | 33 0 53
bdi5 | 42 58 | 33 0 47
bdi8 | 48 52 | 24 0 40
----------------------------------------------------------------------------
欠損数は，Obs=.の部分を確認します。追跡月が進むごとに
欠損が増えていくのが分かります。
次に，欠測パターンをみてみます。
misstable patterns, freq
Missing-value patterns
(1 means complete)
| Pattern
Frequency | 1 2 3 4
-----------------+-------------
52 | 1 1 1 1
|
24 | 1 0 0 0
15 | 1 1 0 0
6 | 1 1 1 0
3 | 0 0 0 0
-----------------+-------------
100 |
Variables are (1) bdi2 (2) bdi3 (3) bdi5 (4) bdi8
この結果から，すべての変数で欠損がない者，つまり，
Patternの欄が全て1である者は，52名しかいないということが
分かります。100名のデータなので，約半数の者で何らかの欠損
があるという状況です。
一方，上記すべての時点が欠損である者は，いちばん下の行の
3名になります（Pattern欄がすべて0）。他にはたとえば，2か月目だけ
データがあり，その後すべて欠損という者が24名（Pattern欄が1,0,0,0），
というように読みます。
この結果から，欠損の代入を何もせずに古典的な反復測定ANOVA
やMANOVAを適用する場合，すべての時点でデータがある，n=52として
解析することになります。 縦断の混合モデル，またはマルチレベル分析等を行う時に，
wide形式のデータセットをlong形式にする必要があります。
適度に欠損値があり，線型混合モデルの解析に適した例
として，RのパッケージHSAUR2のデータセットBtheBを
使ってやってみます。Stataのバージョンは13.1です。
Rからの読み込み方は↓参照
EZRでパッケージ中のRデータセットをStataデータセットに変換する
まずこのままでは，IDがないため，IDを振ります。
egen id=seq()
そして，Stata特有のルールとして，longにする時の
変数名が，言葉+数字という形式になっている必要があるので，
longにまとめるべき抑うつ得点の継時データとなっている変数
の名前を次のように変えます。
rename ( bdi_2m bdi_3m bdi_5m bdi_8m )( bdi2 bdi3 bdi5 bdi8 )
つまり，_とmを取り去った形です。
ここで最初の10名を並べると次のような感じです。
list in 1/10
id bdi2 bdi3 bdi5 bdi8 drug length treatm~t bdi_pre
------------------------------------------------------------
1. 1 2 2 . . No >6m TAU 29
2. 2 16 24 17 20 Yes >6m BtheB 32
3. 3 20 . . . Yes <6m TAU 25
4. 4 17 16 10 9 No >6m BtheB 21
5. 5 23 . . . Yes >6m BtheB 26
------------------------------------------------------------
6. 6 0 0 0 0 Yes <6m BtheB 7
7. 7 7 7 3 7 Yes <6m TAU 17
8. 8 20 21 19 13 No >6m TAU 20
9. 9 13 14 20 11 Yes <6m BtheB 18
10. 10 5 5 8 12 Yes >6m BtheB 20
+------------------------------------------------------------+
そして変換は次のコマンドです。
reshape long bdi, i(id) j(month)
新しく作られる変数が，bdi と month です。
最初の10行をまた出してみます。
list in 1/10
+----------------------------------------------------
id month bdi drug length treatm~t bdi_pre
-----------------------------------------------------
1. 1 2 2 No >6m TAU 29
2. 1 3 2 No >6m TAU 29
3. 1 5 . No >6m TAU 29
4. 1 8 . No >6m TAU 29
5. 2 2 16 Yes >6m BtheB 32
-----------------------------------------------------
6. 2 3 24 Yes >6m BtheB 32
7. 2 5 17 Yes >6m BtheB 32
8. 2 8 20 Yes >6m BtheB 32
9. 3 2 20 Yes <6m TAU 25
10. 3 3 . Yes <6m TAU 25
+----------------------------------------------------
同じidの個人が縦に4行現れるようになりました。また，bdi2からbdi8の変数
はなくなっています。変数bdi2からbdi8の頭の文字であるbdiという変数が新たに
作られ，それぞれの値が順に並んでいます。そして対応する月数の2から8の
部分がmonthという新しい変数の元に縦に並んでいます。他の変数は元の値の
繰り返しです。
メニューからは
Data > Create or change data > Other variable-transformation commands
> Convert data between wide and long
> Convert data between wide and long
と進み，
Long format from wide を選び，
ID variable(s)の欄の右端▼をクリックしてidを選択します。
Subobservation identifier - the j() optionのvariableに
month と入れます。そして
Base (stub) names of X_ij variablesに
bdi と入れてOKです。入力場所は下図参照
この変換したままの状態なら，次のコマンドで元のwide形式に戻す
ことができます。
reshape wide
たくさん描いたヒストグラムを1枚の図にまとめる
必要性があったので，そのやり方を説明します。
おそらく他のグラフでもやり方は同様です。
グラフのまとめ方については以下のコマンドでみられます。
help graph combine
サンプルデータのautoを使って単純な例をやってみます。
sysuse auto
priceとmpgのヒストグラフを別個に出し，それぞれにオプションで
p, mという名前をつけます（何でもよいです）。名前はname( )の中に入ります。
これは「 メモリ中に名前を付けて保存 」というオプションです。
ファイルとして保存することも可能ですが，メモリ中の方が便利です。
hist price, name(p)
hist mpg, name(m)
hist mpg, name(m)
現在メモリ中にあるグラフを一覧にするには
graph dir
とうつと，
m p
と表示されて確認できます。これらを１つにまとめるには
graph combine p m
メモリ中にあるグラフをすべて消したい時は
graph drop _all
です。 たとえば，マルチレベル分析などで，グループごとに平均値や割合を計算して
グループ全員にその平均値の値（例：クラスごとの算数のテストの平均点，
会社の部署ごとのストレス要因得点平均値など），または割合の値（例：
地域ごとの喫煙者の割合，会社の部署ごとの女性比率など）を変数として
追加したい場合があります。
autoデータセットを例に示します。
sysuse auto
で読み込みます。
(1)グループごとの平均値
まず2つのグループだけですが，foreignごとのmpgの平均値を
確認してみましょう。
by foreign, sort: sum mpg
-> foreign = Domestic
Variable | Obs Mean Std. Dev. Min Max
-------------+---------------------------------------------------
mpg | 52 19.82692 4.743297 12 34
------------------------
-> foreign = Foreign
Variable | Obs Mean Std. Dev. Min Max
-------------+---------------------------------------------------
mpg | 22 24.77273 6.611187 14 41
上記結果のように，domesticなら平均19.82692，Foreignなら平均24.77273
となるような新たな変数を作りたいという状況です。次のようにegenを使います。
by foreign, sort: egen mpg2=mean(mpg)
上記コマンドでmpg2という変数で作成できました。確認してみます。
2つのグループを両方表示させられるようにまずpriceで並び替えます。
sort price
list make foreign mpg mpg2 in 1/10
list make foreign mpg mpg2 in 1/10
+--------------------------------------------+
make foreign mpg mpg2
--------------------------------------------
1. Merc. Zephyr Domestic 20 19.82692
2. Chev. Chevette Domestic 29 19.82692
3. Chev. Monza Domestic 24 19.82692
4. Toyota Corolla Foreign 31 24.77273
5. Subaru Foreign 35 24.77273
--------------------------------------------
6. AMC Spirit Domestic 22 19.82692
7. Merc. Bobcat Domestic 22 19.82692
8. Renault Le Car Foreign 26 24.77273
9. Chev. Nova Domestic 19 19.82692
10. Dodge Colt Domestic 30 19.82692
+--------------------------------------------+
foreignのグループ別の平均値が変数として作成されています。
(2)グループごとの割合
単純な例で説明すると，あるグループ内に10名いて，半数が1の値を
とった場合，平均値は (1+1+1+1+1+0+0+0+0+0)/10=0.5なので，割合と
一致します。1の値が3名でも，平均値は (1+1+1+0+0+0+0+0+0+0)/10=0.3
なので割合そのものです。クロス表を描くとすぐに分かります。
tab rep78 foreign, row
Repair |
Record | Car type
1978 | Domestic Foreign | Total
-----------+------------------------------+----------
1 | 2 0 | 2
| 100.00 0.00 | 100.00
-----------+--------------------------------+----------
2 | 8 0 | 8
| 100.00 0.00 | 100.00
-----------+--------------------------------+----------
3 | 27 3 | 30
| 90.00 10.00 | 100.00
-----------+--------------------------------+----------
4 | 9 9 | 18
| 50.00 50.00 | 100.00
-----------+--------------------------------+----------
5 | 2 9 | 11
| 18.18 81.82 | 100.00
-----------+---------------------------------+----------
Total | 48 21 | 69
| 69.57 30.43 | 100.00
つまり，rep78が3だと10%，4だと50%，5だと81.8%，1と2は0%
という結果です。この%を変数として作成するには次のようにします。
詳しくは，Stataの公式FAQsにある通りです。
How can I create variables containing percent summaries?
これも(1)でやったこととほとんど変わらなくて，
最初のbyの後にグループの変数がきて，mean（ ）の所に0-1の2値がくる
というだけです。
by rep78, sort: egen foreignp = mean(foreign)
確認してみます。
sort make
list make rep78 foreignp in 1/10
list make rep78 foreignp in 1/10
+----------------------------------+
make rep78 foreignp
----------------------------------
1. AMC Concord 3 .1
2. AMC Pacer 3 .1
3. AMC Spirit . .2
4. Audi 5000 5 .8181818
5. Audi Fox 3 .1
----------------------------------
6. BMW 320i 4 .5
7. Buick Century 3 .1
8. Buick Electra 4 .5
9. Buick LeSabre 3 .1
10. Buick Opel . .2
+----------------------------------+
rep78には欠損値があるので，ここは無視して他の
foreignpをみると，先ほどのクロス表の結果と一致しています。
このままでもよいですが，出した割合に100をかけてやると，
%と同じ値ができます。
replace foreignp = 100 * foreignp
+----------------------------------+
make rep78 foreignp
----------------------------------
1. AMC Concord 3 10
2. AMC Pacer 3 10
3. AMC Spirit . 20
4. Audi 5000 5 81.81818
5. Audi Fox 3 10
----------------------------------
6. BMW 320i 4 50
7. Buick Century 3 10
8. Buick Electra 4 50
9. Buick LeSabre 3 10
10. Buick Opel . 20
+----------------------------------+
help twoway
でStataでの散布図の描き方解説が見られます。
ここではただの散布図からいろいろ追加情報をいれたものの基本を
例示してみます。
ポイントを一言で言うと，書きたいものを基本のコマンドに
追加していくというだけです。
ちょっと特殊な例ですが，Helpがそうだったので，そのまま説明すると，
共通の横軸になるxの変数，縦軸で異なるy1，y2と
いう変数があった場合，正式な書き方のコマンドの形は，
graph twoway (scatter y1 x) (scatter y2 x)
ですが，省略が色々可能で以下のどれでも同じものが描けます。
twoway (scatter y1 x) (scatter y2 x)
twoway scatter y1 x || scatter y2 x
scatter y1 x || scatter y2 x
scatter y1 y2 x
twoway scatter y1 x || scatter y2 x
scatter y1 x || scatter y2 x
scatter y1 y2 x
以下，3行目の scatter y1 x || scatter y2 x の形式で
進めてみます。まずサンプルデータのautoを呼び出します。
sysuse auto
autoデータの場合，共通のスケールになりそうなy1, y2がなさそう
なので，上のコマンド構造で紹介したものとと少し違いますが，
散布図と回帰直線（lfit）を重ねた図という枠組みでやってみます。
・散布図と回帰直線
scatter mpg weight || lfit mpg weight
・foreignで層別した散布図と回帰直線
scatter mpg weight || lfit mpg weight, by(foreign)
・foreignで層別した散布図と回帰直線+全体のグラフ
scatter mpg weight || lfit mpg weight, by(foreign, total)
1つのグラフにまとめたい時は，ifを使って別個に出していきます。
scatter mpg weight if foreign==0 || scatter mpg weight if foreign==1
回帰直線も入れると長いので，途中に///を入れて改行してます。
以下を実行する際はエディタに貼りつけて実行しないと動きません。
scatter mpg weight if foreign==0 || lfit mpg weight if foreign==0 ///
|| scatter mpg weight if foreign==1 || lfit mpg weight if foreign==1
|| scatter mpg weight if foreign==1 || lfit mpg weight if foreign==1
ピアソンの積率相関では，連続量と連続量の関係ですが，
点双列相関係数(point-biserial correlation)は
2値変数と連続量の関係の相関係数になるみたいです。
とはいっても，出てくる結果は基本的に同じみたい。
Stata13では，2つの平均値の効果量の所から算出できるようになっています。
ためしにやってみます。
sysuse auto
このデータセットでは，foreignが2値データになっています。
まずはピアソンの積率相関から
pwcorr foreign price mpg weight
| foreign price mpg weight
-------------+------------------------------------
foreign | 1.0000
price | 0.0487 1.0000
mpg | 0.3934 -0.4686 1.0000
weight | -0.5928 0.5386 -0.8072 1.0000
次に点双列相関係数をそれぞれ出してみます。
esize twosample price, by(foreign) pbcorr
---------------------------------------------------------
Effect Size | Estimate [95% Conf. Interval]
--------------------+------------------------------------
Point-Biserial r | -.0487195 -.2693882 .1795464
---------------------------------------------------------
esize twosample mpg, by(foreign) pbcorr
---------------------------------------------------------
Effect Size | Estimate [95% Conf. Interval]
--------------------+------------------------------------
Point-Biserial r | -.3933974 -.555367 -.1821459
---------------------------------------------------------
esize twosample weight, by(foreign) pbcorr
---------------------------------------------------------
Effect Size | Estimate [95% Conf. Interval]
--------------------+------------------------------------
Point-Biserial r | .5928299 .4281699 .7051208
---------------------------------------------------------
符号が逆になっていますが，数値は一致しています。 2013/12/29 こまごま修正
Stata13ではη 2 （イータ二乗）とω 2 （オメガ二乗）の効果量（エフェクトサイズ）と
その95%信頼区間が算出可能になっています。
詳しくは，Stata blogの記事
Measures of effect size in Stata 13 を参照
当ブログの記事
1要因分散分析のオムニバス効果量
でRでのやり方を解説したのに習って，Stataでirisデータを読み込んで
解説します。
irisデータのStataでの読み込み方はこちらを参照
Stataでirisデータの散布図を描く
irisデータのcsvファイルを用意して，作業ディレクトリに置いたら，
import delimited iris.csv
またはメニューから
File > Import >Text data (delimited, *csv...)
で読み込み可能。
まず，準備として，品種名の変数speciesは文字データなので，数値に
変換する作業が必要です。文字データを数値データに変換するコマンドは
encode species, gen(nspecies)
で，nspeciesという数値データの変数に変換しました。
次にpetallengthをアウトカムとし，品種間で比べる一要因ANOVAを行うと，
anova petallength nspecies
Number of obs = 150 R-squared = 0.9414
Root MSE = .430334 Adj R-squared = 0.9406
Source | Partial SS df MS F Prob > F
-----------+----------------------------------------------------
Model | 437.102798 2 218.551399 1180.16 0.0000
|
nspecies | 437.102798 2 218.551399 1180.16 0.0000
|
Residual | 27.2225992 147 .18518775
-----------+----------------------------------------------------
Total | 464.325397 149 3.11627783
となります。実は，Stataのhelpの説明によれば，ここで出ているR二乗
R-squared = 0.9414
Adj R-squared = 0.9406
がそれぞれη 2 とω 2 にあたります。デフォルトで情報は出てたみたいですね。
改めてこれらを，信頼区間とともに算出するには上記のanovaコマンドを
実行した後に以下のコマンドを使います。
estat esize
Effect sizes for linear models
-------------------------------------------------------------------
Source | Eta-Squared df [95% Conf. Interval]
--------------------+----------------------------------------------
Model | .9413717 2 .9239033 .9518341
|
nspecies | .9413717 2 .9239033 .9518341
-------------------------------------------------------------------
ω 2 の方は
estat esize, omega
Effect sizes for linear models
-------------------------------------------------------------------
Source | Omega-Squared df [95% Conf. Interval]
--------------------+----------------------------------------------
Model | .9405741 2 .9228679 .9511788
|
nspecies | .9405741 2 .9228679 .9511788
-------------------------------------------------------------------
ためしに線形回帰でもやってみます
reg petallength i.nspecies
Source | SS df MS Number of obs = 150
-------------+------------------------------ F( 2, 147) = 1180.16
Model | 437.102798 2 218.551399 Prob > F = 0.0000
Residual | 27.2225992 147 .18518775 R-squared = 0.9414
-------------+------------------------------ Adj R-squared = 0.9406
Total | 464.325397 149 3.11627783 Root MSE = .43033
------------------------------------------------------------------------------
petallength | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
nspecies |
versicolor | 2.798 .0860669 32.51 0.000 2.627912 2.968088
virginica | 4.09 .0860669 47.52 0.000 3.919912 4.260088
|
_cons | 1.462 .0608585 24.02 0.000 1.341729 1.582271
------------------------------------------------------------------------------
この後にanovaの所で行った
estat esize
estat esize, omega
estat esize, omega
を行っても同じ結果が出てきます。 ※2013/10/26 語句修正 ; 2013/11/07 グループごとの回帰追加 ; 2013/12/23 目次追加，(1)に
factor variablesの確認法，(2)にcontrastコマンド追記 ;
記事の目次
(1)カテゴリ変数のreferenceも結果に出力する方法
(2)1つのカテゴリ変数全体の検定を行ってカイ二乗値とp値を出す
(3)グループごとの回帰（層別）
(2)1つのカテゴリ変数全体の検定を行ってカイ二乗値とp値を出す
(3)グループごとの回帰（層別）
必要だけどつい忘れがちになるポイントのメモです。バージョンはStata 13.0，
使用データセットは，マニュアルのlogisticでexapmle
として使っているHosmer & Lemeshow dataです。インターネットにつながった状態で，
webuse lbw
と打つと読み込まれます。
アウトカムは低体重出生かどうかの2値データ
low
1:birthweight<2500g
0:birthweight>=2500g
説明変数は人種で3区分のカテゴリデータ
race
1:white
2:black
3:other
2:black
3:other
として単純にロジスティック回帰分析を実行すると
logistic low i.race
となり，以下のように
low Odds Ratio Std. Err. z P>z [95% Conf. Interval]
race
black 2.327536 1.078613 1.82 0.068 .9385072 5.772385
other 1.889234 .6571342 1.83 0.067 .9554577 3.735597
_cons .3150685 .0753382 -4.83 0.000 .1971825 .503433
といった結果が出てきます。
i. というのを説明変数の頭につけることで，その変数（ここではrace）を
カテゴリ変数（factor variables）として扱ってくれます。
自動的に複数の2区分変数（ダミー変数）を作成して投入してるイメージです。
デフォルトでは最も低い値がリファレンスになります。
（追記：2013/12/23）
factor variablesがどのように作られているか確認
list race i.race in 1/5
+-----------------------------------+
| 1b. 2. 3. |
| race race race race |
|-------------------------------------|
1. | black 0 1 0 |
2. | other 0 0 1 |
3. | white 0 0 0 |
4. | white 0 0 0 |
5. | white 0 0 0 |
+-----------------------------------+
変数として作成されてはいませんが，このように内部的に
2区分変数（ダミー変数）が作られていることが分かります。
(1)カテゴリ変数のreferenceも結果に出力する方法
デフォルトではオッズ比を出す際にリファレンスとなるカテゴリは出力
されないのですが，オプションにallbaselevelsをつけると，
logistic low i.race, allbaselevels
low Odds Ratio Std. Err. z P>z [95% Conf. Interval]
race
white 1 (base)
black 2.327536 1.078613 1.82 0.068 .9385072 5.772385
other 1.889234 .6571342 1.83 0.067 .9554577 3.735597
_cons .3150685 .0753382 -4.83 0.000 .1971825 .503433
となり，隠れていたwhiteがでてきました。当然のごとく，
オッズ比は1となっています。
また，ラベルをとってどの値と対応しているか確認する場合はnofvlabelを追加します
logistic low i.race, allbaselevels nofvlabel
low Odds Ratio Std. Err. z P>z [95% Conf. Interval]
race
1 1 (base)
2 2.327536 1.078613 1.82 0.068 .9385072 5.772385
3 1.889234 .6571342 1.83 0.067 .9554577 3.735597
_cons .3150685 .0753382 -4.83 0.000 .1971825 .503433
つまり，値1に対する値2,3のオッズ比が算出されているということがわかります。
リファレンスを変えたい場合は，i.raceの項を
ib2.race
ib(last).race
にするとそれぞれ変わります。デフォルトはib(first).raceという意味になっています。
ib(last).race
詳しくは
help factor variables
(2)1つのカテゴリ変数全体の検定を行ってカイ二乗値とp値を出す
それぞれの2区分変数（ダミー変数）の統計量zが出てきますが，論文には
カテゴリ変数全体としてのカイ二乗値を書いてあるのをよく見かけます。
これを出すには，logisticコマンドを実行した後に，さらにtestコマンドを使って，
下記のように結果でパラメータが示されている2区分変数の値を変数の頭につけて，
logistic low i.race
test 2.race 3.race
とします。そうすると，下のような結果が出てきます。
( 1) [low]2.race = 0
( 2) [low]3.race = 0
chi2( 2) = 4.92
Prob > chi2 = 0.0853
つまり，カテゴリ変数raceの全体のワルド検定の結果は
χ 2 =4.92, df=2, p=0.0853
という情報が分かります。
（追記：2013/12/23）
testコマンドよりもっと簡単に全体の検定結果を表示できる
ものがありました。
logistic low i.race
contrast race
contrast race
Contrasts of marginal linear predictions
Margins : asbalanced
------------------------------------------------
| df chi2 P>chi2
-------------+----------------------------------
race | 2 4.92 0.0853
------------------------------------------------
(3)グループごとの回帰（層別）
男女別など，グループ別にロジスティック回帰を行いたい時の方法です。
同じデータを例に，今度は喫煙が低体重出生におよぼす影響を
人種別に解析してみます。まず基本の解析は以下の通りです。
logistic low i.smoke
low Odds Ratio Std. Err. z P>z [95% Conf. Interval]
smoke
smoker 2.021944 .6462989 2.20 0.028 1.08066 3.783112
_cons .3372093 .0724103 -5.06 0.000 .2213694 .5136667
人種別に同じ解析をするには，logisticコマンドの前に以下のような
コマンドを付け加えます。raceの部分を他のカテゴリ変数に変えれば
その変数のカテゴリごとに結果が出ます
by race, sort :
logistic low i.smoke
-> race = white
low Odds Ratio Std. Err. z P>z [95% Conf. Interval]
smoke
smoker 5.757576 3.444621 2.93 0.003 1.782321 18.59916
_cons .1 .0524404 -4.39 0.000 .0357788 .2794949
2013/12/25 こまごま修正と追記
-> race = black
low Odds Ratio Std. Err. z P>z [95% Conf. Interval]
smoke
smoker 3.3 2.775878 1.42 0.156 .6346062 17.16025
_cons .4545455 .2451636 -1.46 0.144 .1579332 1.308222
-> race = other
low Odds Ratio Std. Err. z P>z [95% Conf. Interval]
smoke
smoker 1.25 .8114691 0.34 0.731 .350212 4.461584
_cons .5714286 .1601748 -2.00 0.046 .3298869 .9898259
記事の目次
(1)分散分析表いらない
(2)標準偏回帰係数を出したい
(3)説明変数にカテゴリ変数を使いたい
(4)一連の変数を入れ替えて一気に解析したい
(5)グループごとの回帰（層別）
(2)標準偏回帰係数を出したい
(3)説明変数にカテゴリ変数を使いたい
(4)一連の変数を入れ替えて一気に解析したい
(5)グループごとの回帰（層別）
Stataのサンプルデータセットautoを使って解説します。バージョンはStata 13.1
sysuse auto
まず変数は以下のものを使います。
・走行可能距離の連続量
mpg
・車の重量
weight
・国産か外国産かの2値データ
foreign
0:国産
1:外国産
変数mpgをアウトカムとした重回帰分析は以下のようなコマンドです。
reg mpg weight foreign
実行すると以下のような出力が出ます。
Source | SS df MS Number of obs = 74
----------+---------------------------------- F( 2, 71) = 69.75
Model | 1619.2877 2 809.643849 Prob > F = 0.0000
Residual | 824.171761 71 11.608053 R-squared = 0.6627
----------+------------------------------ Adj R-squared = 0.6532
Total | 2443.45946 73 33.4720474 Root MSE = 3.4071
---------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
----------+----------------------------------------------------------------
weight | -.0065879 .0006371 -10.34 0.000 -.0078583 -.0053175
foreign | -1.650029 1.075994 -1.53 0.130 -3.7955 .4954422
_cons | 41.6797 2.165547 19.25 0.000 37.36172 45.99768
---------------------------------------------------------------------------
(1)分散分析表いらない
メインの回帰係数の部分だけ一覧にして確認したい時は，分散分析表
などあまり関心のない部分は出力しない方がよい時があります。そんな場合は，
reg mpg weight foreign
, noheader
このnoheaderオプションをつけると，上の分散分析表部分は出力されず，
回帰係数の表だけになります。略してnoheadだけでも同様です。
(2)標準偏回帰係数を出したい
オプションにbetaをつけるだけです。
reg mpg weight foreign , beta
また，説明変数とアウトカムすべてを標準化したもので
回帰分析をしても同様の係数が出てきます。
egen zm=std( mpg )
egen zw=std( weight )
egen zf=std( foreign )
reg zm zw zf
egen zw=std( weight )
egen zf=std( foreign )
reg zm zw zf
(3)説明変数にカテゴリ変数を使いたい
年齢層や教育歴，婚姻状態などのカテゴリ変数を使いたい場合は，
factor variablesとなるのでロジスティック回帰分析の際と同じ
です。こちらを参照→ ロジスティック回帰のメモ
reg mpg weight i.foreign , allbaselevels noheader
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight | -.0065879 .0006371 -10.34 0.000 -.0078583 -.0053175
|
foreign |
Domestic| 0 (base)
Foreign | -1.650029 1.075994 -1.53 0.130 -3.7955 .4954422
|
_cons | 41.6797 2.165547 19.25 0.000 37.36172 45.99768
------------------------------------------------------------------------------
(4)一連の変数を入れ替えて一気に解析したい
この方法は回帰分析に関わらず，すべての場面で使えます。
詳しくは→ たくさんの変数に同じ処理を繰り返したいとき 参照
たとえば，下のように説明変数だけを入れ替えて
一連の変数の回帰分析を行いたい時は
reg mpg
weight
, noheader
reg mpg foreign , noheader
reg mpg price , noheader
reg mpg foreign , noheader
reg mpg price , noheader
としますが，数が多くなると書くのが大変になります。
そんな時以下のように書くと，xに代入すべき変数のリストを
varlistの後に書くだけで，同じことができるようになります。
foreach x of varlist weight foreign price {
reg mpg `x', noheader
}
reg mpg `x', noheader
}
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight | -.0060087 .0005179 -11.60 0.000 -.0070411 -.0049763
_cons | 39.44028 1.614003 24.44 0.000 36.22283 42.65774
------------------------------------------------------------------------------
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
foreign | 4.945804 1.362162 3.63 0.001 2.230384 7.661225
_cons | 19.82692 .7427186 26.70 0.000 18.34634 21.30751
------------------------------------------------------------------------------
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
price | -.0009192 .0002042 -4.50 0.000 -.0013263 -.0005121
_cons | 26.96417 1.393952 19.34 0.000 24.18538 29.74297
------------------------------------------------------------------------------
このようにweight, foreign, priceとすべての変数
について入れ替えて実行してくれました。
var1-var100のような書き方もOKです。
(5)グループごとの回帰（層別）
層となる変数にrep78という変数を使用してみます。
以下のようなカテゴリ変数です。
tab rep78,mis
Repair |
Record 1978 | Freq. Percent Cum.
------------+-----------------------------------
1 | 2 2.70 2.70
2 | 8 10.81 13.51
3 | 30 40.54 54.05
4 | 18 24.32 78.38
5 | 11 14.86 93.24
. | 5 6.76 100.00
------------+-----------------------------------
Total | 74 100.00
たとえば，rep78=3の場合のみで解析したい場合，
reg mpg weight if rep78==3, noheader
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight | -.0044096 .0006225 -7.08 0.000 -.0056847 -.0031346
_cons | 33.98076 2.104537 16.15 0.000 29.66981 38.2917
となりますが，
層別にすべての水準で回帰を行いたい場合は以下のようにすると簡単です。
iを1から5の数字に変えて順番に実行してくれます。
forvalues i=1/5{
reg mpg weight if rep78==`i', noheader
}
reg mpg weight if rep78==`i', noheader
}
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight | -.0081081 . . . . .
_cons | 46.13514 . . . . .
------------------------------------------------------------------------------
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight | -.0080464 .0010219 -7.87 0.000 -.010547 -.0055459
_cons | 46.11072 3.45372 13.35 0.000 37.65977 54.56167
------------------------------------------------------------------------------
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight | -.0044096 .0006225 -7.08 0.000 -.0056847 -.0031346
_cons | 33.98076 2.104537 16.15 0.000 29.66981 38.2917
------------------------------------------------------------------------------
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight | -.0049572 .0005596 -8.86 0.000 -.0061435 -.0037709
_cons | 35.89382 1.680188 21.36 0.000 32.33198 39.45566
------------------------------------------------------------------------------
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight | -.01816 .0036908 -4.92 0.001 -.0265091 -.009811
_cons | 69.54448 8.69354 8.00 0.000 49.87832 89.21063
------------------------------------------------------------------------------
（追記： 2013/12/25）
そういえば層別はこちらの方が簡単で分かりやすかったので追記
by rep78, sort: reg mpg weight, nohead
-> rep78 = 1
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight | -.0081081 . . . . .
_cons | 46.13514 . . . . .
------------------------------------------------------------------------------
-> rep78 = 2
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight | -.0080464 .0010219 -7.87 0.000 -.010547 -.0055459
_cons | 46.11072 3.45372 13.35 0.000 37.65977 54.56167
------------------------------------------------------------------------------
-> rep78 = 3
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight | -.0044096 .0006225 -7.08 0.000 -.0056847 -.0031346
_cons | 33.98076 2.104537 16.15 0.000 29.66981 38.2917
------------------------------------------------------------------------------
-> rep78 = 4
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight | -.0049572 .0005596 -8.86 0.000 -.0061435 -.0037709
_cons | 35.89382 1.680188 21.36 0.000 32.33198 39.45566
------------------------------------------------------------------------------
-> rep78 = 5
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight | -.01816 .0036908 -4.92 0.001 -.0265091 -.009811
_cons | 69.54448 8.69354 8.00 0.000 49.87832 89.21063
------------------------------------------------------------------------------
-> rep78 = .
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight | -.0084928 .0022215 -3.82 0.032 -.0155626 -.0014231
_cons | 45.57057 6.414351 7.10 0.006 25.15725 65.9839
------------------------------------------------------------------------------ Stataはバージョン13から標準化効果量（エフェクトサイズ）の算出
に対応しました。2つの平均値の比較については，
CohenのdとHedgesのgが両方とも信頼区間とともに
計算できます。
サンプルデータのautoを使ってコマンド例を示します。
sysuse auto
まず国産車 vs 外国車（foreign）で走行可能距離（mpg）の平均値を比べるt検定は
以下になります
ttest mpg, by(foreign)
Two-sample t test with equal variances
--------------------------------------------------------------------------------------------------------------
Group | Obs Mean Std. Err. Std. Dev. [95% Conf. Interval]
---------+---------------------------------------------------------------------------------------------------
Domestic | 52 19.82692 .657777 4.743297 18.50638 21.14747
Foreign | 22 24.77273 1.40951 6.611187 21.84149 27.70396
---------+-----------------------------------------------------------------------------------------------------
combined | 74 21.2973 .6725511 5.785503 19.9569 22.63769
---------+-----------------------------------------------------------------------------------------------------
diff | -4.945804 1.362162 -7.661225 -2.230384
---------------------------------------------------------------------------------------------------------------
diff = mean(Domestic) - mean(Foreign) t = -3.6308
Ho: diff = 0 degrees of freedom = 72
Ha: diff < 0 Ha: diff != 0 Ha: diff > 0
Pr(T < t) = 0.0003 Pr(|T| > |t|) = 0.0005 Pr(T > t) = 0.9997
次はエフェクトサイズです。ttestの部分を以下のように変えるだけで
簡単に出てきます。
esize twosample mpg, by(foreign)
Effect size based on mean comparison
Obs per group:
Domestic = 52
Foreign = 22
-----------------------------------------------------------------------------
Effect Size | Estimate [95% Conf. Interval]
--------------------+--------------------------------------------------------
Cohen's d | -.9234449 -1.441225 -.3997744
Hedges's g | -.9137865 -1.426151 -.3955932
-----------------------------------------------------------------------------
もちろん，図表などからの要約統計量からでも計算できて，上記の
人数，平均，標準偏差をそれぞれ以下のように入れてもでてきます。
esizei 人数1 平均1 標準偏差1 人数2 平均2 標準偏差2
esizei 52 19.82692 4.743297 22 24.77273 6.611187
Effect size based on mean comparison
Obs per group:
Group 1 = 52
Group 2 = 22
---------------------------------------------------------
Effect Size | Estimate [95% Conf. Interval]
--------------------+------------------------------------
Cohen's d | -.923446 -1.441226 -.3997755
Hedges's g | -.9137876 -1.426153 -.3955942
---------------------------------------------------------
こちらも参照： 平均値，標準偏差，人数の値だけからt検定を行う 自分が調べた範囲ではStata13でもまだこれができないっぽいです。
解決策として，fsumというユーザー作成コマンドをインストールすると
できました。インターネットにつながっている状態で次のようにコマンドを
入力するとインストールできます。
ssc install fsum
ラベルが付いてるサンプルデータで試してみましょう。
sysuse auto
まず，度数分布表では普通にラベルで表示されます。
ここでは，mpgという変数がMileage(mpg)というラベルで
表示されています。
tab mpg
Mileage
(mpg) Freq. Percent Cum.
12 2 2.70 2.70
14 6 8.11 10.81
・・・（略）・・・
35 2 2.70 98.65
41 1 1.35 100.00
Total 74 100.00
次に，summarizeですが，これでは変数名のmpgのみしか出てきません
sum mpg
Variable Obs Mean Std. Dev. Min Max
mpg 74 21.2973 5.785503 12 41
この場合でもラベルを表示させたいとき，fsumに，オプションでlabelをつけると，
以下のように最終列に表示されます。
fsum mpg,label
Variable N Mean SD Min Max日付データは，普通は統計ソフトに読みこんだときは文字型の
mpg 74 21.30 5.79 12.00 41.00 Mileage (mpg)
データになっていると思います*。これを使って日付間の差（年齢算出
や生存日数等）を算出するなどするために，日付変数に変換してStataに
日付のデータだと認識させる必要があります。単純な例を下に
示しました。以下，バージョンはStata12.1での説明です。
//データ作成
set obs 1
gen t1 = "1959/12/1" in 1 //1959年12月1日
gen t2 = "1961/1/1" in 1 //1961年 1月 1日
list t1 t2
set obs 1
gen t1 = "1959/12/1" in 1 //1959年12月1日
gen t2 = "1961/1/1" in 1 //1961年 1月 1日
list t1 t2
+----------------------+
t1 t2
----------------------
1. 1959/12/1 1961/1/1
+----------------------+
この段階ではまだ文字型データです。date関数を使うことで
日付変数に変換できます。
//日付変数に変換 年=Y, 月=M, 日=D
gen t1r = date(t1, "YMD")
gen t2r = date(t2, "YMD")
list t1r t2r
gen t1r = date(t1, "YMD")
gen t2r = date(t2, "YMD")
list t1r t2r
+-----------+
t1r t2r
-----------
1. -31 366
+-----------+
Stataは 1960年1月1日 を起点にしているので，この日を0としてt1が1か月前，
t2が1年後として計算したそれぞれの日付までの日数の数字が上記のように
表示されます。次に，2つの日付間の差を計算します。
gen dif=t2r-t1r
list dif
list dif
+-----+
dif
-----
1. 397
+-----+
これが基本的な日付データの作成方法です。
これから年齢等を計算したい場合は，手動でやるか，
ユーザー作成のpersonageというのを使うとできるみたいです。
findit personage
ところで，実用上日付データには多様な形式があります。
そんな時は，date関数の中身の部分を変えてあげれば対応
できます。
今のStataウィンドウでそのまま続ける場合，
まず前のデータを消去するために
clear
として，新しくデータを入れ直します。
//月-日-年の形式の場合
set obs 1
gen t1 = "12-01-1959" in 1 //12月1日1959年
gen t2 = "01-01-1961" in 1 // 1月 1日1961年
list t1 t2
set obs 1
gen t1 = "12-01-1959" in 1 //12月1日1959年
gen t2 = "01-01-1961" in 1 // 1月 1日1961年
list t1 t2
+-------------------------+
t1 t2
-------------------------
1. 12-01-1959 01-01-1961
+-------------------------+
gen t1r = date(t1, "MDY")
gen t2r = date(t2, "MDY")
list t1r t2r
gen t2r = date(t2, "MDY")
list t1r t2r
+-----------+
t1r t2r
-----------
1. -31 366
+-----------+
と，日付データの並びが違っても，MDYとしてその順番を指定することで，
最初にやった結果と一致しました。月の名前がアルファベットである形式でも
対応します（例：01jan2013 → "DMY"）。
これだけを確認するのなら以下のように打つと簡単です。
display date("01jan2013", "DMY")
19359
*Excelファイルで日付型データにしておいて，StataからExcelファイルを
インポートした場合は自動的に日付変数になってたりします SPSSのデータセット(savファイル)をStataで読み込みたい
場合があるとします。SPSSがあればStataで読める形式に
変換（エクスポート）ができますが，SPSSがない場合はどう
したらよいでしょう？
現段階では，Stata（version 12.1）の方ではsavファイルを
直接インポートする機能はないようです。したがって，少しまわり道
ですが，他の無料で使えるRなどのソフトを使ってデータセットを変換
してからStataで読むという方法が考えられます。
一番直接的なのは， EZR on R Commander を使う方法だと思います。
EZRでは，SPSSデータセットを読めますし，Stataデータセット（dtaファイル）に
直接エクスポートできるからです。しかも，メニューからマウス操作だけで
実行できるので簡単です。
以下，EZR(v1.11)を使って説明します。
下図のように，メニューから
[ファイル]→[データのインポート]→[SPSSのデータセットをインポート]
を選びます。
次に，下図のようなウィンドウが開くので，好きなデータセット名をつけます。
そしてOKを押します。
データセット：＜アクティブデータセットなし＞
→データセット：＜（今名付けたデータセット名）＞
と変わっているでしょうか？
ちゃんと読みこまれたかどうか，[表示]ボタンを押して確認しましょう。
次に，下図のように
[アクティブデータセット]→[アクティブデータセットの更新・保存]
→[アクティブデータセットをエクスポートする（Stata形式）]
→[アクティブデータセットをエクスポートする（Stata形式）]
と選びます。
これで，Stataデータセット（dtaファイル）の完成です。
追記（2013/05/27）：
原因は分かりませんが元のSPSSデータセットが開けなくなったとの報告があったので，
念のため元データのコピーをつくってから試す等，使用は自己責任でお願いします。 変数は英字の方が扱いやすいですが，結果出力時に直観的に把握するために
日本語で名前をつけて，カテゴリ変数のとる値なども日本語に割り当てたいことがあります。
この記事では，変数とその値にラベルをつける方法を説明します。
まず，男性（=0）2名，女性（=1）3名の，計5人のデータを用意してみます。
以下の3行で作成できます。
set obs 5
gen female = 0 in 1/2
replace female = 1 in 3/5
gen female = 0 in 1/2
replace female = 1 in 3/5
データの中身をみてみます
list female
+--------+
female
--------
1. 0
2. 0
3. 1
4. 1
5. 1
+--------+
度数分布表をみると
tab female
female Freq. Percent Cum.
0 2 40.00 40.00
1 3 60.00 100.00
Total 5 100.00
では，まず変数femaleに「性別」というラベルを張ります
label variable female "性別"
tab female
性別 Freq. Percent Cum.
0 2 40.00 40.00
1 3 60.00 100.00
Total 5 100.00
このように，変数名が「性別」にかわりました
次に，値の0,1をそれぞれ男性，女性と置き換えます。Stataは，この作業の時に
それぞれの値に対応するルールを作成して，それを適用するという2つの作業が必要です。
ここでは，ルールの名前はgenderと命名しています。
label define gender 0 "男性" 1 "女性"
次に今作った値のルールgenderを変数femaleに適用します
label values female gender
こうすると，
tab female
性別 Freq. Percent Cum.
男性 2 40.00 40.00
女性 3 60.00 100.00
Total 5 100.00
となっています。元データの方も
list female
+--------+
female
--------
1. 男性
2. 男性
3. 女性
4. 女性
5. 女性
+--------+
とラベルの値が表示されます。
元の数字が見たい時は次のようにオプションをつけます
tab female
, nolabel
性別 Freq. Percent Cum.
0 2 40.00 40.00
1 3 60.00 100.00
Total 5 100.00
となります
【応用】
値のラベルがたくさんあるときは，1行にまとめると混乱するので以下のように書きます
addというオプションをつけて複数行に分けます
label define agecate 0 "20代"
label define agecate 1 "30代" , add
label define agecate 2 "40代" , add
label define agecate 1 "30代" , add
label define agecate 2 "40代" , add
あと値の修正などは，Variables Managerからマウス操作でやった方がやりやすそうな
気がします。
【まとめ】
label variable female "性別"
//変数のラベル作成
label define gender 0 "男性" 1 "女性" //変数の値のラベルルール作成
label values female gender //ラベルルール適用
データハンドリングの中で頻繁に行う作業の備忘録です。
label define gender 0 "男性" 1 "女性" //変数の値のラベルルール作成
label values female gender //ラベルルール適用
Stata付属の有名サンプルデータautoを使ってすぐにでも
試せるように書きました。Stata12.1を使っています。
【連続量をカテゴリに変換，カットオフ分割】
・各カテゴリで人数が等しい割合になるように区分（例：中央値折半，低・中・高の3区分）
・好みの数値で区分（例：抑うつ尺度のカットオフ得点以上）
・好みの数値で区分（例：抑うつ尺度のカットオフ得点以上）
変数作成（拡張版）のegenを使います。詳しくは，マニュアルの
[D] egen - Extensions to generate (p167-)
メニューでは
Data > Create or change data > Create new variable (extended)
からになります。
以下がコマンドです。確認手法を含めて少し冗長に書きましたが，変数作成は
基本的にはegenの行だけです。
sysuse auto
//サンプルデータautoの読み込み
hist price //変数price分布確認
egen cprice = cut(price), group(3)
//等人数になるように低・中・高の3区分でカテゴリ変数cprice作成
tab cprice //作成した変数の内訳確認
egen cprice2 = cut(price), at(0,10000,16000)
//0～10000未満，10000以上16000未満に区分，つまり10000をカットオフ
tab cprice2 //作成した変数の内訳確認
//念のため，cpriceがちゃんと作られてるかどうかの確認
sort price
list price cprice cprice2
hist price //変数price分布確認
egen cprice = cut(price), group(3)
//等人数になるように低・中・高の3区分でカテゴリ変数cprice作成
tab cprice //作成した変数の内訳確認
egen cprice2 = cut(price), at(0,10000,16000)
//0～10000未満，10000以上16000未満に区分，つまり10000をカットオフ
tab cprice2 //作成した変数の内訳確認
//念のため，cpriceがちゃんと作られてるかどうかの確認
sort price
list price cprice cprice2
詳しくは，以下のUCLAのサイトの説明が丁寧です。
How can I recode continuous variables into groups?
【ダミー変数作成(indicator variables)】
ダミー変数とは，ある変数の値の数だけ，それを表す2値変数（0,1）を
作成するものです（ただし回帰分析に使うときはリファレンスにする値の
2値変数を除く）。これは，Stataでは簡単にできます。度数を確かめる
tabコマンドのオプションにgen(newvar)をつけるだけです。
メニューでは
Data > Create or change data > Other variable-creation commands > Create indicator variables
からになります。
以下がコマンドです。ここではcpで始まる変数が3つ作成されます。
tab cprice, gen(cp) //上記で作成したcpriceのダミー変数作成
list cprice cp1 cp2 cp3 //ちゃんとできてるか作成した変数のデータセット表示して確認 まず，下記リンク先，UCLAのStatistical Computingに詳しい解説が
あるので，これに習えば，有名なhsbデータセットで簡単に基本が学べます。
Stataで散布図
http://www.ats.ucla.edu/stat/stata/modules/graph8/twoway_scatter_combine/default.htm
これまでせっかくRのirisデータと親しんだので，上記の方法を
RのirisデータをStataに読ませて再現しようと思います。
irsデータの説明は こちら
まず，Rのirisデータをcsvファイルにエクスポートします。
これは以下の2つのプログラムだけでできます。
setwd(choose.dir())
#これでデータを作成する先を選択します
write.table(iris, "iris.csv", sep=",", row.names=FALSE)
#すでにR内にあるirisデータからiris.csvファイルを作成します。
#最初に設定したフォルダの中にcsvファイルができています。
この方法を使えば，Rのexampleデータセットをなんでもエクスポートして
他のソフトで解析できるので便利です。
ここからはStataです。まずStataを起動します。
次に，csvファイルをインポートします。
メニューからだと
File > Import > Text data created by a spreadsheet
で，読み込むcsvファイルをBrowseから指定した後，
comma-delimited dataにチェックを入れてOKします。
---------------
（追記：2013/12/26）
File > Import >Text data (delimited, *csv...) （Stata13の場合）
---------------
次に，Rでやったように，長い変数の名前を変更します。
メニューからだと
Data > Data utilities > Rename variables
ですが，コマンドだともっと簡単で，（ ）に順に並べるだけです。
なぜかRで表示されていた変数名と違いますが，Stataでは大文字
小文字の区別が（多分）ないのと，変数名でピリオドが認識され
ない，という違いかと思います。
rename (petallength petalwidth species) (pl pw sp)
これで作図準備完了です。
単純には，
twoway scatter pw pl
だけで散布図は描けます。
（1）層別散布図
まずは層別散布図です。カテゴリ変数が品種のspしかないので，これを使います。
twoway (scatter pw pl), by(sp)
なお，byの方にカテゴリ変数を by(female ses)のように追加できる
みたいですが，irisはspしかないのであきらめます。
byオプションを使うと，このように別れて出てきますが，以下のようにすると
Rでやったものと同じように描けました。
（///で複数行に区切ってあるので， Do-file Editorに貼りつけて選択実行 します。）
doedit とコマンドを打つか， アイコン
をクリックしてください。
twoway (scatter pw pl if sp =="setosa") ///
(scatter pw pl if sp =="versicolor" , mcolor(blue)) ///
(scatter pw pl if sp =="virginica" , mcolor(red))
(scatter pw pl if sp =="versicolor" , mcolor(blue)) ///
(scatter pw pl if sp =="virginica" , mcolor(red))
基本的には，scatterを繰り返す感じですね。このままだとラベルがpwの変数ラベル
そのままなので，これを変更するオプションを以下のようにつけます。
twoway (scatter pw pl if sp =="setosa", legend(label(1 "setosa"))) ///
(scatter pw pl if sp =="versicolor" , mcolor(blue) legend(label(2 "versicolor"))) ///
(scatter pw pl if sp =="virginica" , mcolor(red) legend(label(3 "virginica")))
そうすると完成です。
（2）回帰直線
次は線を追加します。
twoway (scatter pw pl) (lfit pw pl)
（3）回帰直線+id表示
まずid変数がないので作成します。これは簡単で，
egen id=seq()
だけでできます。
twoway (scatter pw pl, mlabel(id)) (lfit pw pl)
上のプログラムの最後に , by(sp) とオプションつければ
層別散布図で出るのですが，ごちゃごちゃして見にくかったので
削除しました。 ※2012/04/15 修正
単純な例として，次のようなデータがあるとします。
+---------+
v1 v2
---------
1. 1 0
2. 1 1
3. 1 .
+---------+
v1 v2
---------
1. 1 0
2. 1 1
3. 1 .
+---------+
v1とv2を足し合わせて得点を作りたい場合，どれかの変数に
欠損があっても無視して全部のケースについて平均値を出し，
変数として加えたい時はどうしたらよいでしょう？
次の3通りのコマンドをためしてみます。
gen total1 = v1 + v2 //単純に合計 ---(1)
egen total2 =rowtotal( v1 - v2 ) //関数で合計 ---(2)
egen average =rowmean( v1 - v2 ) //関数で平均 ---(3)
list v1-average
+-------------------------------------+
v1 v2 total1 total2 average
-------------------------------------
1. 1 0 1 1 .5
2. 1 1 2 2 1
3. 1 . . 1 1
+-------------------------------------+
v1 v2 total1 total2 average
-------------------------------------
1. 1 0 1 1 .5
2. 1 1 2 2 1
3. 1 . . 1 1
+-------------------------------------+
(1)と(2)は，まず合計得点を出してみるやり方です。この後変数の数で除す
計算をしてあげれば平均値になります。
(1)は，欠損があるケースは平均値も欠損にする場合有用です
(2)は，欠損値を0に読み変えます。尺度の合計計算には不向きかも知れません。
(3)は，欠損のある変数を抜かして平均値を計算します。
したがって，欠損は関係なく平均値を出したい時は，(3)が有用です。
(3)を使って欠損があったら計算しないようにすることもできます。
egen average2 = rowmean( v1-v2 ) if !missing( v1-v2 )
!missingの中の変数が離れている場合は，(v1,v3)のように , を入れます いつでも参照できるように一覧にしてみます。
適宜，修正・追加していきます。
2012/11/16 追加・修正
--------------------------------------
メモ：
v1など：変数をいれる部分
newvar：新しく作られる変数
--------------------------------------
//ログをとる（解析結果をlogファイルに出力）。ファイルは作業ディレクトリにできる。
//アイコン押すよりこちらの方が楽
log using test
//ここではtestがlogファイルの名前
//この間にプログラム
log close
//この間にプログラム
log close
//データセットを置き換え
use data1.dta, clear
//データセット上書き保存
save data2.dta, replace
//解析の間だけ一時的にプログラム実行（欠損ケースを一時的に抜かして解析したい等）
preserve
//この間にプログラム
restore
//（例：ロジスティック回帰を男女別に）
//男性（sex:0=男性，1=女性，の場合）
preserve
drop if sex==1
logistic incident2 ib(last).agec4
restore
//女性
preserve
drop if sex==0
logistic incident2 ib(last).agec4
restore
//この間にプログラム
restore
//（例：ロジスティック回帰を男女別に）
//男性（sex:0=男性，1=女性，の場合）
preserve
drop if sex==1
logistic incident2 ib(last).agec4
restore
//女性
preserve
drop if sex==0
logistic incident2 ib(last).agec4
restore
//欠損のあるケースを削除
drop if sex==. | age==.
//変数の名前を変える
rename age age2
//名前を変えて新しい変数を作る
gen age2=age
//文字型になっている 文字データ を数値に変換する（文字はラベルとして残る）
/*たとえば，データに性別が「男, 女」で入っているのを「0, 1」に直したい時。*/
encode sex, gen(sex2)
list sex2 //確認
list sex2, nolabel //確認
/*まとめてやりたい時*/
for varlist v1-v3: encode X, gen (Xr)
list sex2 //確認
list sex2, nolabel //確認
/*まとめてやりたい時*/
for varlist v1-v3: encode X, gen (Xr)
//文字型になっている 数字データ を数値型に変換
destring b2 b3, generate(b2n b3n)
//※数字データでも，一部が2~3とか入力されてて文字として読み込まれることがあるので
//ある得点基準によって選ばれた新たな変数を作成する
gen newvar=0
replace newvar=1 if v1>=10 | v2>=11 | v3==12
gen newvar=0 if v1>=2 & v1<=5
replace newvar=1 if v3==1
replace newvar=1 if v1>=10 | v2>=11 | v3==12
gen newvar=0 if v1>=2 & v1<=5
replace newvar=1 if v3==1
//欠損値を数字に変換
recode v1 v2 v3 (.=99)
//値を変換して新たな変数作成
recode v1 v2 v3 (1=0) (2=1), gen(v1r v2r v3r)
//値の変換（接頭語で）
recode v1-v3 (1=0) (2=0) (3=1) (4=1), prefix(dicho_)
//値の範囲の省略
recode v1 (1/2=1) (3/4=2) (5/10=3), gen(v1r)
//たくさんの質問項目の値を一気に変換して新たな変数作成
//文字型変数のたくさんの値から数値型の変数作成（診断名から診断変数作る時とか）
gen dep=0
replace dep=1 if diag=="F320" | diag=="F321" | diag=="F322" | diag=="F323" | diag=="F328" | diag=="F329"
replace dep=1 if diag=="F320" | diag=="F321" | diag=="F322" | diag=="F323" | diag=="F328" | diag=="F329"
//合計得点(欠損は0として変数の合計点)
egen newvar=rowtotal(v1-v3)
//連続変数をカテゴリ分けする（任意の値で）
egen agec=cut(age), at(0,10, 30, 40, 50, 60, 99)
または
gen agec=1 if age~=.
replace agec=2 if age>=30 & age<=39
replace agec=3 if age>=40 & age<=49
replace agec=4 if age>=50 & age<=59
または
gen agec=1 if age~=.
replace agec=2 if age>=30 & age<=39
replace agec=3 if age>=40 & age<=49
replace agec=4 if age>=50 & age<=59
//連続変数をカテゴリ分けする（分割数指定）
egen newvar=cut(v1), group(3)
//ケース1から10までの選択された変数のデータを表示
list id sex age in 1/10
//↑をExcelに出力したいとき （コピーして貼り付けだとうまくいかないので）
export excel id sex age using "test", firstrow(variables) replace in 1/10
//testがファイル名
//グループごとの平均値を作成してデータセットに追加
egen newvar=mean(v1), by(group)
//通し番号(id)を作成
egen id=seq()
sort group
by group: egen ingroupid=seq() //施設ごとに通し番号を付与
追記：2012/10/23
sort group
by group: egen ingroupid=seq() //施設ごとに通し番号を付与
対応のない2群の平均値（mean），標準偏差（sd），人数（n）の情報だけ
からt検定(ttest)を行いたい時があります。これがStataで簡単にできるみたいです。
以下，Stata 12.0で説明します。
メニューから
Statistics > Summaries, tables, and tests > Classical tests of hypotheses
> Two-sample mean-comparison calculator
と進むと，
というウィンドウが出てくるので，これに該当情報を入力する
だけでt検定の結果を出してくれます。とても便利。
コマンドを書くと，もっと簡単にできます
ttesti #obs1 #mean1 #sd1 #obs2 #mean2 #sd2
に，それぞれ情報をいれて実行するだけです。たとえば，
人数1：49，平均1：6.5，標準偏差1：1.2
人数2：51，平均2：7.0，標準偏差2：1.0
というデータがあったら
ttesti 49 6.5 1.2 51 7.0 1.0
とコマンドを打てば，結果が出てきます。参照する結果は以下
t = -2.2671
degrees of freedom = 98
Pr(|T| > |t|) = 0.0256
Rを使う場合はこちら↓
Rを使って平均値と標準偏差と人数だけから2群のt検定（おまけに効果量） Stataで残念な点を挙げるとしたら，出力能力が
あまり多彩でない点でしょうか。自分が勉強足りないだけかもしれませんが，
results windowやlogに出力される方法以外を現時点では知りません。
以下，Stata 12.0で説明します。
出力結果をExcelなどに貼りつけて加工したいときは，
結果を選択して右クリック→Copy Tableでコピーしてから
貼り付けすると大体うまくいきますが，たくさんの出力をいっぺんに
貼りつけたいときには不向きです。
そんなときに使えそうなコマンドが，Ian Watson氏が作成した
コマンドtaboutです。
http://www.ianwatson.com.au/stata.html
Stataでインストールするには，ネットに
つながっている状態で，
findit tabout
と打ち，インストールを選択するとインストールできます。
インストールがすんだら
help tabout
でヘルプが見られます。ここであげられてる例で詳しく説明していきます。
【基本のコマンド例】
exampleのいちばん最初に
. sysuse nlsw88.dta, clear
tabout south race smsa coll [iw=wt] ///
using table2.txt, ///
c(freq row col) f(0c 1p 1p) clab(_ _ _) ///
layout(rb) h3(nil)
とあるのですが，これは複雑でなんだかよくわからなかったので，これを元に
変数やオプションを変えて分かりやすくしてみたのが次のコマンドです。
. sysuse nlsw88.dta, clear
tabout age race ///
using table1.txt, ///
c(freq row ) f(0c 1p) clab(_ _) ///
layout(rb) h3(nil)
ageやraceといったなじみやすそうな変数を選んで，自分の出してみたい
内容に変えてみました。順を追って説明していきます。
試してみたい方は，とりあえずこれらを全部選んで
実行すると，カレントディレクトリにtable1.txtというファイルが作られます
ので，それを開いて，選択してExcelなどに貼りつければきれいな表が
出てきます。結果はresults windowにも出力されます。
カレントディレクトリの位置を知るには，Stata画面の左下に出ていますし，
cd
と打てばパスがresults windowに表示されます。
【コマンド例の解説】
それでは1行ずつのコマンドの意味の説明に入ります。まず，
. sysuse nlsw88.dta, clear
で，Stataに元々入ってるサンプルデータを呼び出します。次に，
tabout age race ///
で，これが本体のコマンドです。ここに，クロス表を出したい変数，
たとえばage raceを入れます。
元々のStataのtabulateで対応してない3要因はもちろんできないです。
///というのは，何行にもわたって1つのコマンドを書いていきたい場合
に必要な記号です。なので，
tabout age race ///
using table1.txt, ///
c(freq row ) f(0c 1p) clab(_ _) ///
layout(rb) h3(nil)
の部分はいっぺんに実行する必要があります。次に，
using table1.txt, ///
ここは，出力するファイル名です。カレントディレクトリ（作業ディレクトリ）
にtable1.txtというファイルが作られます。これより先は" , "で区切られて
いるので，すべてオプションになります。
c(freq row ) f(0c 1p) clab(_ _ _) ///
c(freq row)という所が，出力する内容にあたり，ここでは
度数と行の%を出すという指定にしています。自分はrowを見たい
ことが多いので。縦の%を見たい場合は，colを入れます。
f(0c 1p)というのは，小数点の位置や%の表示を意味していて，
度数は決まって0c，1pは小数点第1位まで表示して，%記号を付ける
という意味です。例えば，f(0c 1p) だと出力が 1,291 9.2%となるという具合。
%をつけたくなければ，f(0c 1)としたら，1,291 9.2となります。
clab(_ _ _)というのは，出力表の3行目の名前を何にするかです。
通常，1行目は変数名，2行目は変数の水準名となります。
3行目には，それぞれの列が何を表わしているかの情報
が入ります。たとえば，No. とか%とかです。ここで使われている _
というのは，空白という意味なんだと思います。さらに，最終行
にある，h3(nil)というのは，そもそも3行目はなしで，という
意味みたいです。
layout(rb) h3(nil)
layout(rb)というのは，表の度数部分と%部分をどういう配置にするかと
いうオプションです。rbだと，度数の表，%の表が縦に出てきます。
layout(cb)にすると，横に連結されます。
h3(nil)というのは上で触れたように，3行目の記述自体をなくす
という指定のようです。ためしに，h3(nil)を削って実行すると，
3行目には空白の行ができてきます。
【繰り返してたくさんの表を出力したい】
たぶんtaboutを使うほどたくさんの処理をする場合には，for varlistとか
で繰り返し処理をする場合がほとんどだと思います。
ここの例で，繰り返してageとその他の変数を変えながらクロス表を
書きたいときには，以下のように書くことでできます。
for varlist race occupation c_city married: tabout age X ///
using table3.txt, ///
append c(freq row ) f(0c 1p) clab(_ _) ///
layout(rb) h3(nil)
for varlist の後に，ageとかけあわせたい変数（ここでは，
race, occupation, c_city, married）を並べ，代入先を X にします。
そして，重要なポイントは， append オプションを追加することです。
こうすると，作成したtable3に繰り返し処理ごとに表を付け加えて
いってくれます。
【その他】
オプションに
style (csv)
と入れて，出力にtable1.csvとすると，直接csvファイルに
出力してくれます。上の例でやると
tabout age race ///
using table1. csv , ///
c(freq row ) f(0c 1p) clab(_ _) ///
layout(rb) h3(nil) style (csv)
追記：2011/12/14
【検定結果も表示させたい】
statsオプションを使ってchi2 gamma V taub lrchi2
といった統計量を表示できるみたいです。カイ二乗検定の結果を表示
させたい場合は，次のオプションを追加します。
stats(chi2)
したがって，上の例でやると，
tabout age race ///
using table1.txt, ///
c(freq row ) f(0c 1p) clab(_ _) ///
layout(rb) h3(nil) stats(chi2)
となります
追記：2012/12/06
上の全部組み合わせたバージョン
（複数の変数のクロス表とカイ二乗検定の結果をcsvファイルに出力）
for varlist sex agec : tabout X stress2 ///
using table1.csv, ///
append c(freq row ) f(0c 2) clab(_ _) stats(chi2) ///
layout(rb) h3(nil) style (csv)
Stataのバージョンは12.0です。
StataでROC (Receiver Operating Characteristic) 曲線を描きたいと思って色々
ネット上の情報を調べた結果，Stataマニュアルが一番分かりやすかったです。
今までSPSSとかで苦労していたので，附属マニュアルが最も分かりやすい
というのは意外でした。学んだことを書いてみます。
y： 病気の有無などの0-1データ
x： スクリーニング尺度の点数データ
とした2つの変数がある場合，
roctab y x, table graph summary
これだけで，0-1データ×尺度の得点のクロス表，ROC曲線のグラフ，AUC
(Area under curve: 曲線下面積)の値と95%信頼区間が全部出てきます。
さらに，
roctab y x, detail
とすると，異なるカットオフポイントごとの感度，特異度，尤度比が出てきます。
ただし，感度，特異度の95%信頼区間が出てきません。かつてはsunspec
というコマンドで出せたらしいのですが，対応がversion8止まりで，11では
動かない感じでした。
調べてみたら， diagt というコマンドをインストールするとできるみたいです。
上記のxを求めたいカットオフポイントで区切った0-1変数（x2）に直して，
diagt y x2
と打つと，ちゃんと感度，特異度の95%信頼区間が出力されました。
roctabと同様の働きをする他のコマンドでlrocというのもありますが，こちらより
roctabの方が簡単だし，出てくる情報も多いようです。
実際にデータを使って試してみたい場合は，Stataマニュアル
のExamplesにあるコマンドで体験できます。 webuse hanley と打つだけで
データセットが読み込まれます。多分ネットにつながってないとだめだろうけど。
メニューでは下記のように選んでいくとたどりつきます
Statistics > Epidemiology and related > ROC analysis >
Nonparametric ROC analysis without covariates
※2013/05/12 Rの部分は間違えてたので修正中