レジュメ9 DIDの基礎と応用
1 因果推論入門
ヘドニック・モデルを用いて、水害リスクが地価に与える影響を分析します。2015年9月7日に発生した台風第18号の影響により、同年9月10日に茨城県常総市において鬼怒川の堤防が決壊し、大規模な浸水被害が発生しました。この浸水被害により、これまで水害リスクについて正しく認知されていなかった地域において、水害への考え方が大きく変わった可能性があります。そこで、そのような水害リスク認知の変化を、地価の変化を通して検証します。このように、ある原因によって、どのような結果が得られたのか分析することを因果推論と言います。
鬼怒川決壊の様子（出所 国土交通省(2015)）
1.1 欠落変数バイアスの問題
まず、水害リスクが地価に与える影響について検証するために、ある年の地価公示データを用いて、下記のようなシンプルな重回帰分析を考えてみましょう。
\[
\begin{align}
&ln(地価_i)=\alpha_0+\beta_0*浸水想定区域ダミー_i+
\beta_{1i}*土地属性_{1i}+・・・+ \beta_{ki}*土地属性_{ki}+e_i ・・・(1)
\end{align}
\]
\(浸水想定区域ダミー_i\)は浸水想定区域内であれば1、そうでなければ0の値をとるダミー変数、\(土地属性_{1i},・・・,土地属性_{ki}\)は標準地iの土地属性（面積、容積率、駅までの距離など）、\(e_i\)は誤差項です。ここで関心のある係数は、浸水想定区域ダミーの係数である\(β_0\)です。これは、浸水想定区域外の地域に比べて、水害リスクの高い浸水想定区域内では、平均的にどの程度地価が異なるのかを示しています。
しかし、このようなクロスセクションデータによるシンプルな重回帰分析では、欠落変数バイアスの問題が生じる可能性が高く、OLS推定量にはバイアスが含まれてしまいます。例えば、浸水想定区域に指定されている地域は河川に近いため、日当たり・眺望の良さなどがある一方、地盤が弱く、液状化や地盤沈下のリスクが高いといった特徴があります。これらの特徴は地価に大きく影響するため、地域要因を十分にコントロールしない限り、浸水想定区域ダミーと誤差項が相関してしまい、欠落変数バイアスの問題が生じます。
では、あらゆる地域要因を数値化して、説明変数に追加すれば解決できるかというとそれも難しいです。この地域要因には、教育・文化環境、地域のつながりなどの観測不可能、または観測が非常に困難な要因も含まれるため、重回帰分析で地域要因をコントロールするには限界があります。
川沿いの住宅のメリット（出所 アットホーム）
1.2 パネルデータ分析
パネルデータ分析では、各観測個体の時間を通じた変化を捉えることができるため、観測可能、不可能に関わらず、各観測個体の時間を通じて一定な要素を取り除いた推定が可能となります。ここで、欠落変数バイアスの原因となる要素は、時間を通じて一定であると仮定すると、パネルデータ分析を用いることで欠落変数バイアスを解消することができます。
そこで本演習では、地価公示データをパネルデータ化し、パネルデータ分析の基本である固定効果モデルの推定方法を学びます。具体的には、以下のヘドニック・モデルを推定します。
\[
\begin{align}
&ln(地価_i)=\alpha_i+\beta_1*(shinsui_i*after_t)+
\beta_2*(higai_i*after_t)+u_{it} ・・・(2)
\end{align}
\]
\(α_i\)は標準地の固定効果、\(shinsui_i\)は浸水想定区域内かつ浸水被害のない地域であれば1、そうでなければ0の値をとるダミー変数、\(higai_i\)は浸水想定区域内かつ浸水被害のあった地域であれば1、そうでなければ0の値をとるダミー変数、\(after_t\)は被害後ダミーで公表年度が2016年以降であれば1、そうでなければ0の値をとるダミー変数、\(u_{it}\)は誤差項です。なお、定数項は\(α_i\)に含まれます。
この推定モデルは、固定効果モデルと呼ばれます。固定効果モデルとは、「欠落変数バイアスの原因となる変数は時間を通じて一定」という仮定の下、通常の重回帰モデルに固定効果\(α_i\)を追加したモデルです。\(α_i\)は\(i\)ごとに異なる時間を通じて一定な要素をまとめたもので、本演習で使用する地価公示データは観測単位が標準地であるため、標準地レベルで時間を通じて一定な要素となります。したがって、前節で採用していた時間を通じて一定であると考えられる変数（地積、容積率、都心主要駅までの距離など）は、固定効果\(α_i\)として処理しています。また、観測不可能な時間を通じて一定な要素も\(α_i\)に含まれます。
このように\(α_i\)を追加することで、誤差項\(u_{it}\)には、説明変数と相関するであろう時間を通じて一定の要素は含まれないため、誤差項と説明変数は無相関となり、欠落変数バイアスは解消されると考えます。
1.3 差の差推定（difference-in-differences estimation; DID）
その一方で、浸水想定区域ダミー、浸水地域ダミーを定数項ダミーとして採用すると、固定効果\(α_i\)として処理されてしまい、本来関心のある影響をみることができません。そこで、被害後ダミーとの交差項を採用します。これは、浸水被害を自然実験とみなし、浸水想定区域外を対照群、浸水想定区域と浸水地域を介入群とした差の差推定（difference-in-differences
estimation; DID）を行っていることになります。
差の差推定とは、実験データが得られない状況で、因果関係を推定する手法の一つです。ここで、浸水地域（介入群）に属する標準地iの被害前の地価を\(chika_{T,i0}\)とし、被害後の地価を\(chika_{T,i1}\)とします。すると直感的には、\(E(chika_{T,i1}-chika_{T,i0})\)を浸水被害が地価に与える平均的な影響だと考えがちですが、これは誤りです。例えば、この地域では浸水被害前から住環境が年々悪化し、地価が下がり続けているという地域のトレンドがあった場合、被害前後の平均的な地価の差は、浸水被害の影響＋地域のトレンドを捉えてしまいます。そこで、浸水地域に近接する地域のデータを利用します。近接地域（対照群）に属する標準地iの被害前の地価を\(chika_{C,i0}\)とし、被害後の地価を\(chika_{C,i1}\)とし、近接地域における被害前後の平均的な地価の差
\(E(chika_{C,i1}-chika_{C,i0}
)\)を定義します。これは、浸水被害を受けなかった場合の地価の変化、つまり地域のトレンドと考えられます。
ここで、「被害がなければ浸水地域と近接地域との間で地域のトレンドは同じである」という仮定（平行トレンドの仮定）が成り立つとしましょう。すると、\(E(chika_{T,i1}-chika_{T,i0} )-
E(chika_{C,i1}-chika_{C,i0})\)によって、浸水被害の影響のみを捉えることが可能となります。このように、介入群の前後差から対照群の前後差を引くことで介入効果を推定する手法であることから、差の差推定（difference-in-differences
estimation; DID）と呼ばれます。
また、差の差推定は、介入群ダミー\(（T_i）\)、介入後ダミー\(（after_i）\)、介入群ダミー×介入後ダミー\((T_i×after_i)\)を採用した線形回帰モデルによっても推定が可能となります。
\[
\begin{align}
&Y_i=\alpha+\beta_1*T_i+\beta_2*after_t+\beta_3*(T_i*after_t)+u_{it}
・・・(3)
\end{align}
\]
ここで、\(α\)は対照群の介入前の\(Y_i\)の平均値、\(α+β_2\)は対照群の介入後の\(Y_i\)の平均値、\(α+β_1\)は介入群の介入前の\(Y_i\)の平均値、\(α+β_1+β_2+β_3\)は介入群の介入後の\(Y_i\)の平均値を示しています。したがって、対照群の前後差は\(β_2\)、介入群の前後差は\(β_2+β_3\)となり、差の差をとった\(β_3\)が介入効果となります。
1.4 Two-ways fixed effect model（TW-FE）
２式は、浸水被害を自然実験とみなし、対照群を浸水想定区域外、介入群を浸水想定区域内かつ浸水被害のない地域、浸水想定区域内かつ浸水被害のあった地域としたモデルです。このモデルを推定することで、浸水被害が地価に与える影響を推定でき、その影響は、水害リスク認知の変化として解釈することが可能です。もし、浸水想定区域の設定が適切で、事前から水害リスクが十分に認知されていれば、浸水被害があったとしても水害リスクの認知は変わらず、地価に影響はないと考えられます。したがって、\(β_1\)と\(β_2\)を推定することによって、浸水被害前から水害リスクの認知が十分であったかを検証することができます。
\[
\begin{align}
&ln(地価_i)=\alpha_i+\beta_1*(shinsui_i*after_t)+
\beta_2*(higai_i*after_t)+u_{it} ・・・(4)
\end{align}
\]
差の差推定によって因果関係を推定する前提として、平行トレンドの仮定が満たされている必要があります。そこで、モデル２では被害後ダミーの代わりに、2012年度を基準とする公表年度ダミー
\(year_t\)を採用することで、浸水被害前のトレンドが介入群と対照群で異なるのか検証します。また、浸水被害の影響は一時的ではなく、長期にわたって影響し続ける可能性があるため、被害後の長期的な影響も検証します。このように、固定効果\(α_i\)と時点効果\(year_t\)を導入したモデルは、two-ways fixed
effect model（TW-FE）と呼ばれます。
1.5 Rによる推定
1.5.1 Rの設定
Rを用いて、固定効果モデルを推定します。まずRを起動します。使用するパッケージをインストールします。（既にインストールしてある場合は、この作業は必要ありません。）
#install.packages("sf")
#install.packages("car")
#install.packages("tidyverse")
#install.packages("estimatr")
#install.packages("modelsummary")インストールしたパッケージを呼び出します。
#library(sf)
#library(car)
#library(tidyverse)
#library(estimatr)
#library(modelsummary)setwd()を用いて、ワーキングディレクトリを設定します。（getwd()を実行すると、現在のディレクトリを確認できます。）
#setwd("C:/gis_ouyou/flood/") 1.5.2 データの読み込みと抽出
茨城県の県西地区、県南地区の住宅地のデータのみを用います。そこでまず、県西地区、県南地区の名前をまとめたベクトルを作成します。
#city <- c("かすみがうら","つくば","つくばみらい","下妻","五霞","利根","取手","古河","土浦","坂東","守谷","常総","桜川","牛久","石岡","稲敷","筑西","結城","茨城八千代","茨城境","茨城河内","阿見","龍ケ崎") 次に、st_read関数を使用して「chika_ibaraki_final.shp」を読み込み、filter関数を使用して県西地区と県南地区の住宅地のデータのみを抽出し、subset関数とselect関数を使用して必要な変数のみを残したものを「chika」というオブジェクトに保存します。
#chika <- st_read(dsn="chika_ibaraki_final.shp") %>%
#filter(L01_001 == "000",L01_022%in%city) %>%
#subset(select=c(L01_085, L01_086, L01_087, L01_088, L01_089,
#L01_090,L01_091, L01_092, L01_093, L01_094, d_shinsui, d_higai))Point：パイプ演算子（% > %）は、左側の出力やデータを右側の関数の第一引数にするものです。
1.5.3 パネルデータの作成
ここで、データには標準地固有のIDが付与されていないため、mutate関数を使用してIDを付与します。
#chika <- mutate(chika,id = row_number())変数のL01_085からL01_094までが、2012年から2021年の地価（円／㎡）ですが、データを見ると左寄せになっており、文字列として認識されています。そこで、mutate_at関数を用いて、これらの変数を数値に変換します。
#chika <- mutate_at(chika, vars(L01_085:L01_094), as.numeric)地価公示は標準地の正常な価格を示したもので、この標準地は代表性、中庸性、安定性、確定性の４つの観点から選定されています。毎年、４つの点に合致しているか点検されており、不適格な標準地については選定替が行われています。そのため、パネルデータ化するにあたって、2012年から2021年の間で選定替が行われた標準地を除きます。そこで、filter関数を使用して2012年から2021年の間で地価の値が0のデータを除きます。
#chika <- filter(chika, L01_085>0,L01_086>0,L01_087>0,L01_088>0,L01_089>0,L01_090>0,L01_091>0,L01_092>0,L01_093>0,L01_094>0)次に、colnames関数を使用して、各年の地価の変数名を西暦に変更します。
#colnames(chika) <- c("2012","2013","2014","2015","2016","2017","2018","2019","2020","2021","d_shinsui","d_higai","geometry","id")パネルデータのデータ形式は、ワイド形式とロング形式があります。ここまで作業してきたchikaデータはワイド形式になっていますが、パネルデータ分析をするためには、ロング形式にする必要があります。そこで、gather関数を使用してワイド形式からロング形式に変換したchika_longを作成します。
#chika_long <- gather(chika, key=year, value=price, "2012","2013","2014","2015","2016","2017","2018","2019","2020","2021")1.5.4 変数の作成
さらに次のように、新たな変数を作成します。①地価を自然対数に変換します。②ifelse関数を使用して、浸水想定区域ダミー（浸水想定区域内かつ浸水被害のない地域であれば1、そうでなければ0の値をとるダミー変数）を作成します。③被害後ダミー（2016年以降のデータであれば1、そうでなければ0の値をとるダミー変数）を作成します。④浸水想定区域ダミー×被害後ダミー、浸水地域ダミー×被害後ダミーを作成します。ここで、浸水想定区域を越えて浸水した地域はないため、浸水地域ダミーは、実質的に浸水想定区域内かつ浸水被害があった地域であれば1、そうでなければ0となるダミー変数となります。
#chika_long$ln_price <- log(chika_long$price)
#chika_long$d_nshinsui <- ifelse(chika_long$d_shinsui==1 & chika_long$d_higai==0,1,0)
#chika_long$d_after <- (chika_long$year>=2016)
#chika_long$nshinsui_after <- (chika_long$d_after)*(chika_long$d_nshinsui) #chika_long$higai_after <- (chika_long$d_after)*(chika_long$d_higai)年ダミーはforループを使用して、次のように作成します。
#for(i in 2012:2021){ #chika_long[sprintf("d_%d", i)] <- ifelse(chika_long$year==i, 1, 0) #}また、浸水想定区域ダミー×年ダミー、浸水地域ダミー×年ダミーもforループを使用して、次のように作成します。
#for(i in 2012:2021){ #chika_long[sprintf("nshinsui_%d", i)] <- ifelse(chika_long$d_nshinsui==1 & chika_long$year==i, 1, 0)
#chika_long[sprintf("higai_%d", i)] <- ifelse(chika_long$d_higai==1 & chika_long$year==i, 1, 0) #}1.5.5 モデルの推定
これで推定作業の準備が整いました。まず、モデル１をlm_robust関数を使用して推定します。結果を見たい場合はsummary関数を使用します。
#base_fe <- lm_robust(ln_price ~ nshinsui_after+higai_after+d_after,
#data = chika_long,
#clusters = id,
#fixed_effects = ~ id,
#se_type = "stata") #summary(base_fe)続いて、モデル２をlm_robust関数を使用して推定します。
#long_fe <- lm_robust(ln_price ~ nshinsui_2013+nshinsui_2014+nshinsui_2015+nshinsui_2016+nshinsui_2017+nshinsui_2018+nshinsui_2019+nshinsui_2020+nshinsui_2021+higai_2013+higai_2014+higai_2015+higai_2016+higai_2017+higai_2018+higai_2019+higai_2020+higai_2021+d_2013+d_2014+d_2015+d_2016+d_2017+d_2018+d_2019+d_2020+d_2021,
#data = chika_long,
#clusters = id,
#fixed_effects = ~ id,
#se_type = "stata")
#summary(long_fe)最後に推定結果の表を作成します。表をまとめるパッケージはいくつかありますが、lm_rubustについては、前節で紹介したstargazerパッケージが対応していないため、ここではmodelsummaryパッケージを使います。まず、２つ以上のモデルを並べて表記したい場合は、モデルをまとめたリストを作成します。また、変数名を編集したい場合は、変数名のベクトルを作成しておきます。最後に、msummary関数を使用して、推定結果をきれいにまとめます。
#base_model <- list(base_fe,long_fe) #var_names <- c('nshinsui_after' ='浸水想定区域ダミー×被害後ダミー','higai_after' ='浸水地域ダミー×被害後ダミー','nshinsui_2013' ='浸水想定区域ダミー×2013年ダミー','nshinsui_2014' ='浸水想定区域ダミー×2014年ダミー','nshinsui_2015' ='浸水想定区域ダミー×2015年ダミー','nshinsui_2016' ='浸水想定区域ダミー×2016年ダミー','nshinsui_2017' ='浸水想定区域ダミー×2017年ダミー','nshinsui_2018' ='浸水想定区域ダミー×2018年ダミー','nshinsui_2019' ='浸水想定区域ダミー×2019年ダミー','nshinsui_2020' ='浸水想定区域ダミー×2020年ダミー','nshinsui_2021' ='浸水想定区域ダミー×2021年ダミー','higai_2013' ='浸水地域ダミー×2013年ダミー','higai_2014' ='浸水地域ダミー×2014年ダミー','higai_2015' ='浸水地域ダミー×2015年ダミー','higai_2016' ='浸水地域ダミー×2016年ダミー','higai_2017' ='浸水地域ダミー×2017年ダミー','higai_2018' ='浸水地域ダミー×2018年ダミー','higai_2019' ='浸水地域ダミー×2019年ダミー','higai_2020' ='浸水地域ダミー×2020年ダミー','higai_2021' ='浸水地域ダミー×2021年ダミー')
#msummary(base_model,statistic = 'std.error', stars = c("*" = .05, "**" = .01, "***" = .001),coef_map=var_names,title = "推定結果", 'markdown')