gravity.xlsx データには、2005 年の 1年間の世界の2国間の貿易額が収録されている。
Rでデータを読み込む。
library(readxl)
gravity <- read_excel("gravity.xlsx")
head(gravity)## # A tibble: 6 × 20
## impor…¹ expor…² year imports gdp_e…³ gdp_im…⁴ join_…⁵ join_…⁶ expor…⁷ impor…⁸
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 AFG ABW 2005 0 NA 6.81e 9 NA NA 533 4
## 2 AGO ABW 2005 0 NA 3.06e10 NA 1994 533 24
## 3 ALB ABW 2005 0 NA 8.38e 9 NA 2000 533 8
## 4 ANT ABW 2005 4.34e+3 NA NA NA NA 533 530
## 5 ARE ABW 2005 9.51e-1 NA 1.33e11 NA 1994 533 784
## 6 ARG ABW 2005 6.58e-1 NA 1.83e11 NA 1967 533 32
## # … with 10 more variables: contig <dbl>, comlang_off <dbl>, colony <dbl>,
## # dist <dbl>, REPlandlocked <dbl>, PARTlandlocked <dbl>, religion <dbl>,
## # onein <dbl>, bothin <dbl>, nonein <dbl>, and abbreviated variable names
## # ¹importer, ²exporter, ³gdp_exporter, ⁴gdp_importer, ⁵join_exporter,
## # ⁶join_importer, ⁷exporternum, ⁸importernum変数 importer は輸入国の ISO コードと呼ばれる 3桁の国コードである。同様に、exporter は輸出国 の ISO コードである。例えば、日本の ISO コードは JPN、アメリカの ISO コードは USA である。また、ISOコード（exporter, importer）の他に、輸入額（imports）、輸出国のGDP（gdp_exporter）、輸入国のGDP（gdp_importer）、輸出国と輸入国の間の距離（dist）が含まれている。
R で次のような重力方程式を推定することを考える。
\[ 輸入額 = 定数 \times \frac{輸出国のGDP^{\alpha} \times 輸出国のGDP^{\beta}}{輸出国と輸入国の間の距離^{\gamma}} \times e ^{\delta \times 言語の共通性} \] この式は非線形のためコンピューターの能力上、推定が難しい。そのため、従来は両辺の対数をとって、線形にされてきた。線形にされた式の推定は最小二乗法で比較的簡単に推定できる。そこでまず、回帰分析に用いる変数の対数を取る。Rのコードは次の通りである。
gravity$limports<-log(gravity$imports)
gravity$lgdp_exporter<-log(gravity$gdp_exporter)
gravity$lgdp_importer<-log(gravity$gdp_importer)
gravity$ldist<-log(gravity$dist)ここで、0より大きい値しか対数値にできないことに注意が必要である。 そのため、輸入額が0より大きい値のサンプルを作成しておく。
gravity2<-subset(gravity,gravity$imports>0)通常、貿易データは不均一分散の性質を持つので、デフォルトのlmで推定するのではなく、lm_robustで推定する。そのために必要な パッケージ{estimatr}をインストールする。
install.packages("estimatr")パッケージ{estimatr}の使い方は、公式ページで確認できる。
不均一分散頑健な標準誤差を計算するようlm_robustで重力方程式を推定するコードは以下の通りである。summary(ols)により回帰分析の結果が出力される。
library(estimatr)
ols<-lm_robust(limports~lgdp_exporter+lgdp_importer+ldist+comlang_off,data = gravity2)
summary(ols)##
## Call:
## lm_robust(formula = limports ~ lgdp_exporter + lgdp_importer +
## ldist + comlang_off, data = gravity2)
##
## Standard error type: HC2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) -33.751 0.352954 -95.63 0.000e+00 -34.4433 -33.0596 19973
## lgdp_exporter 1.226 0.007648 160.29 0.000e+00 1.2108 1.2408 19973
## lgdp_importer 0.951 0.007571 125.61 0.000e+00 0.9361 0.9658 19973
## ldist -1.374 0.020435 -67.24 0.000e+00 -1.4141 -1.3340 19973
## comlang_off 1.293 0.050422 25.65 8.613e-143 1.1945 1.3922 19973
##
## Multiple R-squared: 0.6424 , Adjusted R-squared: 0.6423
## F-statistic: 9877 on 4 and 19973 DF, p-value: < 2.2e-16なお、回帰分析は、以下のコードでも実行できる。
# ols<-lm_robust(gravity2$limports ~ gravity2$lgdp_exporter + gravity2$lgdp_importer + gravity2$ldist + gravity2$comlang_off)ここで説明されているように、stargazerではなく、Texregやmodelsummaryを用いた方が簡単に、lm_robustを用いた不均一分散頑健な推定結果を表示できる。
#install.packages("modelsummary")
library(modelsummary)
#(1) lmを使用して通常の回帰分析を実行し、ols1に保存。
ols1 <- lm(limports~lgdp_exporter+lgdp_importer+ldist+comlang_off,data = gravity2)
#(2) lm_robustを使用して回帰分析を実行し、不均一分散頑健な標準誤差を計算し、ols2に保存。
ols2<-lm_robust(limports~lgdp_exporter+lgdp_importer+ldist+comlang_off,data = gravity2)
modelsummary(list(ols, ols2))| (1) | (2) | |
|---|---|---|
| (Intercept) | −33.751 | −33.751 |
| (0.353) | (0.353) | |
| lgdp_exporter | 1.226 | 1.226 |
| (0.008) | (0.008) | |
| lgdp_importer | 0.951 | 0.951 |
| (0.008) | (0.008) | |
| ldist | −1.374 | −1.374 |
| (0.020) | (0.020) | |
| comlang_off | 1.293 | 1.293 |
| (0.050) | (0.050) | |
| Num.Obs. | 19978 | 19978 |
| R2 | 0.642 | 0.642 |
| R2 Adj. | 0.642 | 0.642 |
| AIC | 93181.8 | 93181.8 |
| BIC | 93229.2 | 93229.2 |
| RMSE | 2.49 | 2.49 |
より手の込んだ推定結果表も作成できる。 参考：modelsummary: regression tables
modelsummary(
list(ols, ols2),
estimate = c("estimate",
"{estimate}{stars}"),
coef_omit = "Intercept",
gof_omit = 'DF|Deviance|AIC|BIC|RMSE',
notes = list('* p<0.1; ** p<0.05; *** p<0.01')
)| (1) | (2) | |
|---|---|---|
| lgdp_exporter | 1.226 | 1.226*** |
| (0.008) | (0.008) | |
| lgdp_importer | 0.951 | 0.951*** |
| (0.008) | (0.008) | |
| ldist | −1.374 | −1.374*** |
| (0.020) | (0.020) | |
| comlang_off | 1.293 | 1.293*** |
| (0.050) | (0.050) | |
| Num.Obs. | 19978 | 19978 |
| R2 | 0.642 | 0.642 |
| R2 Adj. | 0.642 | 0.642 |
| * p<0.1; ** p<0.05; *** p<0.01 |
なお、stargazerを用いた表の作成は、stargazerがlm_robustに対応していないためそのままではできない。
ここで説明されているような工夫を施すことで、stargazerを用いた表の作成ができるようになる。
まず、lmを使用して通常の回帰分析を実行する。そのあとで、starprep関数を用いて、標準誤差を不均一分散頑健な標準誤差に変換する。
library(stargazer)##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazer#(1) lmを使用して通常の回帰分析を実行し、ols1に保存。
ols1 <- lm(limports~lgdp_exporter+lgdp_importer+ldist+comlang_off,data = gravity2)
#(2)starprep関数を用いて、ols2の標準誤差を不均一分散頑健な標準誤差に変換して、ols2の推定結果を表示する。
stargazer(ols1, type="text",se = starprep(ols2))##
## ================================================
## Dependent variable:
## ----------------------------
## limports
## ------------------------------------------------
## lgdp_exporter 1.226***
## (0.008)
##
## lgdp_importer 0.951***
## (0.008)
##
## ldist -1.374***
## (0.020)
##
## comlang_off 1.293***
## (0.050)
##
## Constant -33.751***
## (0.353)
##
## ------------------------------------------------
## Observations 19,978
## R2 0.642
## Adjusted R2 0.642
## Residual Std. Error 2.492 (df = 19973)
## F Statistic 8,969.064*** (df = 4; 19973)
## ================================================
## Note: *p<0.1; **p<0.05; ***p<0.01#(3)通常のOLSの標準誤差は以下で表示できる。
stargazer(ols1, type="text")##
## ================================================
## Dependent variable:
## ----------------------------
## limports
## ------------------------------------------------
## lgdp_exporter 1.226***
## (0.008)
##
## lgdp_importer 0.951***
## (0.008)
##
## ldist -1.374***
## (0.022)
##
## comlang_off 1.293***
## (0.051)
##
## Constant -33.751***
## (0.350)
##
## ------------------------------------------------
## Observations 19,978
## R2 0.642
## Adjusted R2 0.642
## Residual Std. Error 2.492 (df = 19973)
## F Statistic 8,969.064*** (df = 4; 19973)
## ================================================
## Note: *p<0.1; **p<0.05; ***p<0.01