の続きです。今回は課税所得と昼夜人口比率の関係を調べてみます。前回、前々回は昼夜人口比率と人口密度、昼夜人口比率と商品販売額と1つの変数と1つの変数の関係でしたが、今回は、人口密度と商品販売額をコントロールした上で、課税所得と昼夜人口比率の関係を調べようと思います。
具体的には、人口密度と商品販売額が同じとき、昼間の人口のほうが多い(day == 1)都市と、夜の人口のほうが多い都市(day == 0)では、課税所得に違いはあるかどうか？です。
tax = β0 + β1 * day + β2 * mitsudo + β3 * hanbai + u
という線形モデルを想定して、lm()関数でOLSで係数を推定します。
summary()関数で係数を確認します。
tax = 2257 -237 * day + 0.055 * mitsudo + 0.23 * hanbai + u
という推定式です。人口密度と商品販売額が同じならば、昼間の人口のほうが多い都市は、そうでない都市と比べて23万7千円ほど課税所得が少ない、ということですね。
p値は0.0192と5%よりも小さい水準です。
人口密度、商品販売額ともに係数は正の値なので、どちらも値が大きいほど、課税所得も大きいということです。これは常識的な感覚と一致しますね。
OLSでの推計で気を付けるのは、残差が均一分散かどうかです。グラフで確認してみます。
右上のグラフを見ると、均一分散していないように見えます。
lmtestパッケージのbptest()関数でチェックしてみます。
p値が2.247e-09とかなり小さい値です。ということは残差は不均一分散しているということです。
OLSの推計で残差が不均一分散のときは、ロバスト標準誤差を求めるか、WLS(Weighted Least Square)で係数を推定するか、です。
まずはロバスト標準誤差を求めてみます。
ロバスト標準誤差でもdayのp値は0.02138と5%よりも小さな値ですから、やっぱり昼夜人口比率の大きい都市と小さい都市では課税所得に違いがあるという結論はかわらないですね。
WLSでも係数を推計してみましょう。
Step1でOLSでの残差の2乗の対数を求めます。
Step2でそれを回帰分析して、
Step3で、その回帰分析の推定値から重みを求め
step4で、その重みを考慮して回帰分析します。
Step5で、結果をsummary()関数で表示しています。
WLSでの推計結果もOLSと大きく違いは無いですね。dayの係数のp値がさらに0に近くなりました。
最後に、stargazer()関数で、OLSとWLSの推計結果をきれいな表にします。
OLSでもWLSでも、昼間の人口のほうが多い都市は、23万円ぐらい課税所得が少ないです。
今回は以上です。
はじめから読むには、
です。
今回のコードは以下になります。
#
# taxとdayの関係(mitsudo, hanbaiをコントロール)
# 線形モデル
lm_mod <- lm(tax ~ day + mitsudo + hanbai, data = df)
#
# lm_modのサマリー
summary(lm_mod)
#
# lm_modのプロット
par(mfrow = c(2, 2))
plot(lm_mod)
par(mfrow = c(1, 1))
#
# Heteroskedasticityの有無
library(lmtest)
bptest(lm_mod)
#
# ロバスト標準誤差
library(sandwich)
coeftest(lm_mod, vcov = vcovHC(lm_mod, type = "HC1"))
#
# WLS(Weighted Least Square)
# Step 1 残差の2乗の対数(正の数にするため)
logu2 <- log(residuals(lm_mod)^2)
#
# Step 2 logu2の推定モデル
varreg <- lm(logu2 ~ day + mitsudo + hanbai, data = df)
#
# Step 3 重みを生成
w <- 1 / exp(fitted(varreg))
#
# Step 4 WLSモデル
wls_mod <- lm(tax ~ day + mitsudo + hanbai, data = df,
weights = w)
#
# Stpe 5 wls_modのサマリー
summary(wls_mod)
#
# lm_mod, wls_modの比較
library(stargazer)
stargazer(lm_mod, wls_mod, type = "text")
#
(冒頭の画像は Bing Image Creator で生成しました。プロンプトは Close-up of Japanese beautyberry (Callicarpa japonica), under the blue high sky, photo です。
