Stata

Stataで作業を行うにあたって、ログを取っておくことは非常に重要です。というのも、研究の再現性を確保するためには過去に自分が何を行ったかわかっていないといけないからです。また、結果が長すぎるとStataの画面で見ることが難しいケースもあるため、logを取っておくとメモ帳で見ることができて便利です。

基本的な方法は、ログをとり始めるときにlog using ファイル名、logを閉じるときにlog closeと書きます。

以下にログの取り方の例を記載しておきます。

毎日日付を書き替えるのが面倒なので、2,3行目で今日の日付を表すマクロを作成することで対処しています。

cd "c:/desktop"

local todaysdate = c(current_date)

local todaysdate = subinstr("`todaysdate'", " ", "", .)

log using "projectA\log\ProjectA_log_`todaysdate'.log", replace

....

log close

.dtaファイルなど、ファイルを異なるフォルダへ移動させたいが記録を残したい場合、stata上で行うことができる。

copy 元のファイルのパス 新しく作成するファイルのパス

で可能となる。例えば、downloadフォルダにあるauto.dtaというファイルをdocument/dtaに入れたい場合、下のコードになる

copy "C:/download/auto.dta" "C:/document/dta/auto.dta"

コピーした上で元のファイルを消せば、実質的に移動させたことになる。

erase ファイルのパス　

で消えるので、

erase "C:/download/auto.dta"

とすればよい。

改行の方法

コードが長くなった時に、読みやすさのために改行したいことがよくあります。Stataには主に3通りの改行方法があります。

1.  /*  */を使う

Stataでは、/*以降の文字は、/*まで無視されます。

ですので、reg x1 /*ここはメモ*/ x2 はdoファイルで実行するとreg x1x x2と同じです。

この性質を応用して、改行したい行の末尾に/*を付け、その次の行の初めに*/を付けると改行できます。

個人的にはこれを一番使用しています。

例：

reg x1 x2 /*

*/ x3, r

est sto result1

predict u, residual

2. ///を使う

改行する行の末尾に///を付けるだけです。おそらくこれが一番一般的に使用されています。

reg x1 x2 ///

 x3, r

est sto result1

predict u, residual

3. #delimit; を使う

#delimit;でこのコマンドが開始して#delimit crで終了します。

このコマンドが生きている間は、改行を;で行います。あまり使ったことがありませんが、まれにこのコマンドを使ったほうが見やすい場合もあります。推奨はしません。

#delimit ; 

reg 

x1 

x2 

x3; 

est sto result1 

#delimit cr //ここで#delimitはおわり 

predict u, residual

コメントアウトの方法

コメントアウトとは、コードには書いているが実際にコードを動かしたときに実行されない部分を指します。

上でもコード内に登場していますが、コメントアウトは3つ方法があります。

1. 　行の先頭に*を付ける

これは、１行を丸ごとコメントアウトしたいときに使います。

2. 　行の末尾に//を付ける

行の末尾に簡単なコメントを付けたいときやオプションを消したいときに使えます

3.  /*  */を使う

複数行にわたるコメントや、コード内の一部をコメントアウトしたいときに使えます。

四捨五入したいときは、round(x,y)でxをyで丸めてくれる。例えば、round(15,10)だと20になる。yを書かないと小数第一位で四捨五入する。小数第一位で四捨五入したいときは、int(x)を用いることもできる。

切り捨てはfloor(x)で行うが、どの位で切り上げをするか選択できないので、例えば百の位を切り上げたいときには1000で割ってからfloor()を使ってその後に1000掛けるとよい。

切り上げは、ceil(x)を用いる。

割り算の余りを表示するときは、mod(a,b)を使用します。例えば、mod(2134,100)=34です。

上から通しで番号をつけたいときはgen id=_nもしくはegen id=seq()

n個ずつでつけたいときは、egen id=seq(), b(n)

例えば、egen id=seq(), b(3)だと上からidは111222333と続く

また、aからbまでを繰り返したいときはegen id=seq(), f(a) t(b)

例えば、123123123...としたいときはegen id =seq(), f(1) t(3)

文字列化から数字へ変換したいときに使われるのがdestringである。新しい変数を作りたいときは、destring var, gen(newvar)、元の変数名のままだとdestring var, replaceとなる。

数字から文字列にしたいときはtostringであり、オプションは上のdestringと同じである。

destringでは数字以外のものが入っているとエラーが起こるが、この原因を見つけるときに役に立つのが

gen nonnumeric =real(var)==. である。実数であれば0、そうでなければ1のダミー変数を作成してくれる。

それぞれのグループごとに変数を作成したいとする。例えば、市町村別の人口データがあったときに都道府県の人口を作成したいとする。

変数は都道府県コード(prefcode)、市町村コード(citycode)、年(year)、人口(pop)とする。

方法①sum()関数を使う

bysort prefcode year: egen totalpop =sum(pop)

新しい変数を作るので、元の変数は残っている。

方法②　collapse関数を使う

collapse (平均や和など実行したい演算) 集約したい変数, by(グループ化する基準の変数)

なので、今回は下のようなコードになります。

collapse (sum) pop, by(prefcode year)

もともとの市町村別の人口は消える。また、コードに記載していない変数(今回はcitycode)が消える。

mvencode 変数名, mv(0)で欠損値をゼロに変えることができる。

全ての変数でやりたいときは、

mvencode _all, mv(0)

また、特定の変数のみを変えたいときにfor文を使うことも可能である。

foreach var of varlist gdp* pop* {

replace `var' = 0 if `var' =.
}

ここでのvarというのはローカルマクロであるため、プログラムを走らせている間のみ有効である。

この例ではgdpで始まる変数とpopで始まる変数の欠損値をすべてゼロにしている。

また、gdp*という変数がgdp2000,gdp2001, ... gdp2020などとなっている場合は以下のコードでも可能である。

foreach i of numlist 2000/2020{

replace gdp`i' =0 if gdp`i'=.
}

また、ローカルマクロを使わない方法として

forvalues i=2000/2020{

replace gdp`i' =0 if gdp`i'=.

}

も可能である。

しかし、個人的には２つめのコマンドのほうが好きである。理由は、foreachの後には文字列を、forvaluesの後には数字を書くようにすると覚えやすいからである。

また、数字のところを2000/2020と書いているがこれは2000から2020まで1ずつ数を増やして実行していく。これに対して、偶数年のみに実行したいとすると2000(2)2020と書けばよい。括弧の中が前の数字に足していく数とイメージしてもらえばよいので、2000/2020と2000(1)2020は同じである。

ストリング型の変数についてダミー変数を作成したいとする。

たとえば、prefという変数に都道府県名が入っているとしてdpref1, dpref2, ... dpref47まで47つダミー変数を作成する。（実際に分析で使用するダミー変数は47-1=46個）

この場合は、tab pref, gen(dpref)と入力することで可能となる。

回帰分析などで使用するときはすべて書くのは面倒なので、reg y $xlist dpref* のようにワイルドカード(*)を用いることで簡略化できる。

別の方法で、egen dpref = group(pref)とすると1から47のカテゴリー変数が作成される。その後reg y $xlist i.dpref　とすると上と同様の結果となる。

続いて、数値の場合を考える。

例えば年齢を表す変数ageがあり、20歳未満を1とするyoungというダミー変数を作成したいとする。そして、欠損値については.としたい。このような場合には、

gen young=0

replace young=1 if age<20

replace young=. if age=. // if missing(age)でも可

とする。

また、上で作成したdprefという1から47の都道府県コードが入ったカテゴリー変数について考える。

上のreg y $xlist i.dprefようにダミー変数を作成しないという方法も存在しているが、xi: reg y $xlist i.dprefでダミー変数を自動的に作成してから回帰分析を行ってくれる。この方法だと_Idpref_2から_Idpref_47という名前のダミー変数を作成してもらえる。

これと関連して、cond()関数も紹介する。

cond(条件,a,b) は条件が正しいとa、正しくないとbとなる。

下のようなデータがあり、欠損値に上の値を代入したいとする。

id year var

1 2000 2

1 2001 .

2 2000 5

2 2001 5

3 2000 3

3 2001 .

このような場合に上の値を代入したい、つまり一つ目の欠損値に"2"、二つ目に"3"を入れたいとする。

そのようなケースでは、下のコードを用いる。

replace var=var[_n-1] if missing(var)

[_n-1]が一つ上の変数という意味で、if missing(var)がvarという変数が欠損値の場合という意味である。

細かなことはこのサイトを参照のこと。

inlistとregexmについて

Stataでは"もしくは"を表すときに|を用いるが、これをいくつも連続して使うとしんどい。例えば、prefという都道府県を表す変数があったとして東北地方をdropしたいとする。

このようなケースでは、通常だと

drop if pref =="Aomori"| pref=="Iwate" | pref=="Miyagi" | pref=="Akita" | pref=="Yamagata"| pref=="Fukushima"

というようなコードを書く。

これだとめんどくさいので今回はinlist()関数を導入する。

drop if inlist(pref, "Aomori", "Iwate", "Miyagi", "Akita", "Yamagata", "Fukushima")

これで上と同じ意味になる。

ちなみにinlist()関数単体で使えば、下記のように東北地方なら1のダミー変数の作成に使える。

gen tohoku = inlist(pref, "Aomori", "Iwate", "Miyagi", "Akita", "Yamagata", "Fukushima")

また、かなり特殊な例として山の入っている都道府県を除外したい場合はregexm()関数を用いる。

drop if regexm(pref, "yama") |  regexm(pref, "Yama")と入力すると名前に山の入った都道府県(山形、富山、山梨、和歌山、岡山)がドロップされる。

今回のケースでも、gen yama =regexm(pref, "yama") |  regexm(pref, "Yama")と入力することで上記の県ならば１のダミー変数が作成できる。

if文内での省略方法

また、ダミー変数の作成については上記で述べたようにinlist()やregexm()が使えるほかに、generateで右辺に条件式を持ってくるだけで作成できる。

例えば、価格(price)が100未満のものを1とするダミー変数dpriceを作成する場合

gen dprice =0

replace dprice =1 if price<100 

と書くのがふつうであるが、

gen dprice =price<100

としても同じ結果となる。

一般的な読み込み方は、import excel using "ファイル名", firstrow sheet("sheet名")

firstrowオプションは、エクセルの1行目に変数が入っているとき（通常はそのようなデータが多い）に1行目を変数名として認識するものである。

そして、1行目にタイトルが入っていたりするデータや2列目からデータが入っているような場合は、
読み込む範囲を指定して、import excel using ファイル名, cellrange(B2:K100)などとする。

そして、エクセルファイル内にシートが複数存在している場合は、まず

import exvel using ファイル名, describe と入力することでエクセルシートの範囲やsheet名の一覧を見ることが可能である。

つづいて、return list と入力することでn番目のsheet名がr(worksheet_n)、範囲がr(range_n)というローカルマクロを作成してくれる。

そこで、10番目のsheetを読み込みたいときは、import excel using "ファイル名", sheet(`r(worksheet_10)')とすればよい。

注意点としては、マクロを書く際には`'で囲うということである。

さらに、文字化けする場合はエンコーディングを変更するとよい。オプションでencoding(shift-jis)をつけるとなおることが多い。

結果を表やグラフなどきれいな形で出力したい場合、まずは結果をStata内に保存してその後に作業を行うことが必要である。

そこで、まずは

estimate store 保存したい名前

で分析結果を保存する。

その後に結果を出力するためのコードを書く。

１．分析結果を表で表すとき

esttabコマンドを用いて結果を出力する。

基本的にはesttab 保存した名前 で結果を出力することができるが、もう少し見た目を整えるためにオプションを付ける。

２． 分析結果をグラフで表示するとき

coefplotコマンドを用いる。こちらも、coefplot 保存した名前 で結果を出力できるが、オプションを用いて見た目を整えるのが一般的である。

３．オプションについて

3.1 必要な係数のみ出力したいケース

結果を出力する際に、分析で用いたすべての変数を表やグラフで表すのではなく主要な結果のみを記すのが一般的である。

そのためのオプションとしてkeep()オプションが存在する。

例えば、reg y x_1 x_2 x_3 x_4という分析を行いestimate store basicmodel と結果を保存したとする。そして、主要な結果はx_1とx_2であるためにこれらのみを出力したいとする。

このようなケースにおいては、esttab basicmodel, keep(x_1 x_2)と入力することで必要な表が手に入れられる。

また、drop()オプションを使うことで不要な変数を除くことができる。上記と同じ結果を得ようとすればesttab basicmodel, drop(x_3 x_4 _cons)と入力すればよい。_consは定数項で、Stataにおいては普通に回帰分析をすれば自動で入ってくる。定数項を抜きで分析を行いたい場合は、noconstantオプションを用いることで可能となる。

こんな感じの結果になる。

----------------------------------------------------

                              (1)             (2)   

                            Price           Price   

----------------------------------------------------

x_1                        -63.21          -52.22   

                          (-0.75)         (-0.62)   

x_2                         2.442***        2.111** 

                           (3.55)          (3.41)   

----------------------------------------------------

Observations                   69              69   

----------------------------------------------------

t statistics in parentheses

* p<0.05, ** p<0.01, *** p<0.001

3.2　交差項を選択して出力したいケース

Stataでは、回帰分析の際に#を用いることで交差項を作成することができる。

例えば、連続変数のx_1およびx_2を用いて交差項を作成したい場合はc.を各変数の前に付けて間に#を入れることで作成ができる。つまり、c.x_1#c.x_2とすることでx_1*x_2の交差項が作成できる。さらに#をもう一つ加えることでx_1、x_2、そしてx_1*x_2の3つを表すことができる。

すなわち、reg y x_1 x_2 c.x_1#c.x_2 とreg y c.x_1##c.x_2は同じである。

連続変数の交差項を選択して結果を出力する場合は、*(ワイルドカード)を用いる。

例えば、c.x_1#c.x_2の結果を出力したい場合は、オプションにてkeep(*.x_1#*.x_2)とすればよい。

つづいて、d_1やd_2というダミー変数の交差項を作成したい場合はc.の代わりにi.を用いることで可能となる。

そして、この結果の一部を出力したい場合は、iの代わりに出力したいカテゴリーを選んで

keep(1.d_1#1.d_2)のようにすればよい。

こんな結果になる。

------------------------------------

                              (1)   

                            Price   

------------------------------------

x_1 # x_2                  -0.133   

                          (-1.83)   

------------------------------------

Observations                   69   

------------------------------------

t statistics in parentheses

* p<0.05, ** p<0.01, *** p<0.001

3.3 係数の名称を変えたいとき

labelオプションを使うことで、結果にラベルを表示させることが可能となる。

例：

label var x_1 mpg

label var x_2 weight

esttab model1, keep(x_1 x_2) label

--------------------------------------------

                      (1)             (2)   

                    price           price   

--------------------------------------------

mpg                -63.21          -52.22   

                  (-0.75)         (-0.62)   

weight              2.442***        2.111** 

                   (3.55)          (3.41)   

--------------------------------------------

N                      69              69   

--------------------------------------------

t statistics in parentheses

* p<0.05, ** p<0.01, *** p<0.001

3.4 有意水準の星を変えたいとき

個人的に使っているのはstar(* 0.1 ** 0.05 *** 0.01)オプションである。アスタリスク以外も指定できるし、有意水準の数字も自由に変えられる。

デフォルトは、star(* 0.5 ** 0.01 *** 0.001)となっている。

esttab model1, keep(x_1 x_2) label star(* 0.1 ** 0.05 *** 0.01)

--------------------------------------------

                      (1)             (2)   

                    price           price   

--------------------------------------------

mpg                -63.21          -52.22   

                  (-0.75)         (-0.62)   

weight              2.442***        2.111***

                   (3.55)          (3.41)   

--------------------------------------------

N                      69              69   

--------------------------------------------

t statistics in parentheses

* p<0.1, ** p<0.05, *** p<0.01

3.5 タイトルをつけたいとき

mtitle()オプションで可能となる。例えば、esttab model1 model2, mtitle("NO FEs" "TWFE" )などとすることで、下のような結果となる。

--------------------------------------------

                      (1)             (2)   

                    NO FE            TWFE   

--------------------------------------------

mpg                -63.21          -52.22   

                  (-0.75)         (-0.62)   

weight              2.442***        2.111***

                   (3.55)          (3.41)   

--------------------------------------------

N                      69              69   

--------------------------------------------

t statistics in parentheses

* p<0.1, ** p<0.05, *** p<0.01

3.5 結果をエクスポートしたいとき

using を使う。例えば、esttab model1 model2 model3 using example.csv, keep(x1 x2)

分析結果のグラフとして、例えばMikawa(2025)のFig.5のようなグラフを作成したいとします(この「折りたたみ可能なグループ」へ見やすい大きさでの貼り付け方がわかりませんでした)。

この場合、使用するコマンドはcoefplotです。例えば、よく使用されるauto.dtaを使用した分析例を示します。

たとえば、以下のように回帰分析の結果を係数の点と信頼区間のバーで表すことができます。変数名にはラベルが表示されています。

sysuse auto

qui reg price mpg rep78 weight length

coefplot

ただし、このコマンドではすべての係数が１つの表にまとまってしまっており、やや見づらくなっています。そこで表示する変数を以下のように指定してあげることが可能です。

coefplot, keep(mpg rep78 weight length)

また、切片を消したい場合には

coefplot, drop(_cons)

としてやっても同様の結果です。

さらに、現在は横軸に係数、縦軸に変数名となっていますがverticalオプションを付けることで逆にできます。

coefplot, drop(_cons) vertical 

注意点として、xtregやreghdfeコマンドではcoefplotを使用できません。

基本統計量などの結果を用いて計算を行いたい場合を考える。

よくある基本統計量の出し方としてsummarizeがあるが、この結果は以下のようになっている。

データは、よく使われるautoを使用する。

. sysuse auto

(1978 Automobile Data)

. sum price

    Variable |        Obs        Mean    Std. Dev.       Min        Max

-------------+---------------------------------------------------------

       price |         74    6165.257    2949.496       3291      15906

. display `r(max)'-`r(min)'

12615

上記のように、summarizeで求められた結果は、r(なまえ)という名前のマクロで保存されている。

return listと打てば下のように一覧が返ってくる。詳しくは、こちらを参考にしていただきたい。

. return list

scalars:

                  r(N) =  74

              r(sum_w) =  74

               r(mean) =  6165.256756756757

                r(Var) =  8699525.974268788

                 r(sd) =  2949.495884768919

                r(min) =  3291

                r(max) =  15906

                r(sum) =  456229

結果を表示したいときはdisplayコマンドを使用する。例えば、最大値を表示したいときはdisplay r(max)と入力すれば結果が返ってくる。

この結果を用いて計算等を行う場合に一番最初に思いつくのはgenerateコマンドだが、これはベクトルである。（シートの1列すべて同じ値で追加される）かわりにscalarを用いることでスカラーを作成し、使用するメモリを節約することができる。

scalar max =r(max)とすることでmaxという定数が作成できる。この後で、最大値からの差という変数divを作成したいとすると、generate div = max- priceとすることで可能となる。または、マクロのまま使用してgenerate div = `r(max)'- priceとすることも可能である。

基本的にはreghdfeを使うようにしましょう。ssc install ftools とssc install reghdfe でインストールできます。メリットとしては、主にxtsetしなくてよいこと、決定係数を出してくれること、固定効果のレベルをabsorb()オプションで複数選択できることなどがあります。

xtregを使用した場合、固定効果が個別効果のみになってしまい、時間効果が入っていません。ですので、two-wayの固定効果をとろうとすると、

xtreg y x1 x2 ... i.time, feとする必要があります。

詳しくは慶應の松浦先生のnoteを参照してください。

Sun and Abraham (2020)を行うことができるeventstudyinteractの使い方
https://asjadnaqvi.github.io/DiD/docs/code/06_eventstudyinteract/

その他のStaggered DIDについても、上と同様にAsjad Naqviさんのサイトにあります

また、日本語では日経学会の2023年春季大会で慶應義塾大学の小西先生が発表された資料がとても分かりやすいです。

日付について、たまに持っているデータによっては2022/6/20などのような文字で入っている場合がある。その場合は、そのままではStataは認識してくれないので、date関数を用いる。

例えば、上に書いたような文字列がdayという変数に入っていて日付をあらわすためにtimeという変数を作成したいとする。この場合は、

gen time= date(day, "YMD", .)

とする。意味について説明する。date関数の中で2つのカンマで区切られて3つの単語(?)が入っている。意味は、dayという変数の中にはYear, Month, Dateの順番で日付を表す変数が入っているということである。注意点としては、""を用いないとエラーが出ることである。時間が入っている場合は、h,m,sを用いて表す。例えば、2022/6/20 12:59:48などと入って入れば"YMDhms"と書く。このアルファベットは繰り返しになるが、元の変数に登場する順番で書かないといけない。3つめの単語はカンマだけだと1960年1月1日から何日後かを返す。年が2桁の場合にはカンマ以外の数字を用い、書いた年を超えない最大の年が返される。別の方法として、"MD19Y"だと1900年代が返ってきたり"MD20Y"だと2000年代が返ってきたりする。

例えば、gen time= date("99/9/25", "YMD", 1900)の場合には、1899年9月25日となる。

この時は-22012という数字が返ってくるが、これだけでは使えないので日付の形にしたい。

そこで、

format time %tc

と入力することで25sep1899というように日付になる。

capture:　

capture: を前につけることでもしエラーメッセージが出るような場面でもエラーを出さずに続けることができる。doファイルを使用している場合には便利な機能である。

例えば、新しくログを開始したいがこれまでログを開いていたのか分からないとき次のように書けばよい。

capture: log close

erase     

erase "ファイル名"で不要なファイルを消すことができる

例：eraseold_version.dta

以下工事中

まずはshapeファイルをStataで使用するためにdtaへ変換する方法です。

Stata15以降であればspshape2dtaコマンドを使用することで可能となります。

spshape2dta ファイル名

と入力することでファイル名.dtaとファイル名_shp.dtaという2つのファイルが作成されます。

なお、ポリゴンデータの場合は重心を_CX(経度)_CY(緯度)で出してくれます。

ちなみに、これをすると毎回日本語が文字化けして困っています。

解決策がunicode translate データ名.dtaです。clearしてunicode encoding set shiftjisと書いてから実行しましょう。

例：

cd c:/user/desktop \\ディレクトリ指定

unicode encoding set shiftjis \\shiftjisへエンコーディング変更

clear

unicode translate XXX.dta\\文字化け解消

Stataで２地点間の距離を測定したいときはgeodistを使います。使い方は、geodist 緯度1 経度1 緯度2 経度2, generate(変数名)です。注意点としては、arcGISなどで作業するデータはxが経度でyが緯度のことが多く順番が逆になるということです。

Stataでジオコーディングを行う方法としてgeocode3などがある。

ジオコーディングとは住所から緯度経度を求める方法であるが、日本で対応しているものが少ない。国内だと例えばCSVアドレスマッチングサービスを用いるという方法でCSVファイルを作成することができる。

Getting Started with Stata for Windows

初心者向けのマニュアル

What's it like–Getting started in Stata 

初心者向けのStataのチュートリアルビデオ

ネコでもわかるStata入門

東大の別所さんが昔作ったもの

Stataのグラフについて

The Mixetape

　誤植が多いですが、文字による説明が中心で読みやすいです

　Stataのコードに加えてRとPythonのコード例も掲載されています