1 はじめに

2022年1月・2月号の「実証ビジネス・エコノミクス」の第6回「戦略は企業の特性が決める：静学ゲームの推定［応用編］」に関するRプログラムの紹介である。本プログラムによって誌面で紹介した分析を再現しつつ、実際に手を動かしながら学ぶことを目的とする。

1.1 留意事項

本連載の内容、およびサンプルコード等の資料は、情報の提供のみを目的としていますので、運用につきましては十分にご確認をいただき、お客様ご自身の責任とご判断によって行ってください。これらの情報の運用結果により損害等が生じた場合でも、日本評論社および著者はいかなる責任を負うことはできませんので、ご留意ください。

1.2 全体の流れ

はじめに
Rに関する下準備
データの準備
企業レベルの変数の記述統計
Berry (1992)による推定
推定結果
反実仮想分析

1.3 謝辞

今回のプログラムの作成に際して島本幸典さん(東京大学大学院経済学研究科)にご尽力頂きました。また、プログラムの実行確認に際して、小磯慎士さん(早稲田大学政治経済学部)にご協力頂きました。この場を借りて、御礼申し上げます。

2 Rに関する下準備

以下で読み込んでいるパッケージについて、まだインストールしていない場合には、install.packages()でインストールを行う。

# ワークスペースを初期化
rm(list = ls())

# パッケージを読み込む
require(tidyverse)

## Loading required package: tidyverse

## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --

## v ggplot2 3.3.3     v purrr   0.3.4
## v tibble  3.1.1     v dplyr   1.0.5
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   1.4.0     v forcats 0.5.1

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

require(tictoc)

## Loading required package: tictoc

require(skimr)

## Loading required package: skimr

require(modelsummary)

## Loading required package: modelsummary

require(kableExtra)

## Loading required package: kableExtra

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

require(doParallel)

## Loading required package: doParallel

## Loading required package: foreach

## 
## Attaching package: 'foreach'

## The following objects are masked from 'package:purrr':
## 
##     accumulate, when

## Loading required package: iterators

## Loading required package: parallel

require(doRNG)

## Loading required package: doRNG

## Loading required package: rngtools

3 データの準備

# 病院データを読み込む
data_raw <- readr::read_csv(file = "data_hospital_chapter6.csv", 
                            locale = locale(encoding = "UTF-8"))

## 
## -- Column specification --------------------------------------------------------
## cols(
##   .default = col_double(),
##   Prefecture = col_character(),
##   SecondaryAreaName = col_character(),
##   Management = col_character()
## )
## i<U+00A0>Use `spec()` for the full column specifications.

# ダミー変数を作成
data_cleaned <- data_raw %>% 
  # 大学（公立、国立、私立大学）ダミー
  dplyr::mutate(DaigakuDum = if_else(
    Management %in% c('公立大学法人','国（国立大学法人）','私立学校法人'),1, 0)) %>% 
  # 病床数0ダミー
  dplyr::mutate(ZeroBedDum = if_else(NumBeds == 0, 1, 0))

以下のように単位を調整する。
- NumBeds: 100病床
- Population: 100万人
- Menseki: 100km2
- TaxableIncome: 100万円
また、NumBedsとPopulation、TaxableIncomeは対数変換した変数を作成する。

# 単位の調整と、対数変換した変数の作成
data_cleaned <- data_cleaned %>% 
  dplyr::mutate(
    NumBeds = NumBeds/100,
    LogNumBeds = log(NumBeds+0.01),
    Population = Population/10^6,
    Menseki = Menseki/100,
    TaxableIncome = TaxableIncome/1000,
    LogPop = log(Population),
    LogIncome = log(TaxableIncome))

4 企業レベル変数の記述統計

MRIを購買した病院と購買していない病院の記述統計

data_cleaned %>% 
  dplyr::group_by(MRIOwnDum) %>% 
  dplyr::select(MRIOwnDum, Kyukyu, Kinou, Sien, Hyoka, 
         DepNeurology ,DepNeurosurgery, NumBeds,
         ZeroBedDum, DaigakuDum) %>% 
  skimr::skim() %>% 
  skimr::yank("numeric") %>% 
  dplyr::select(skim_variable, MRIOwnDum, mean, sd) %>% 
  as.data.frame() %>% 
  tidyr::pivot_wider(names_from = MRIOwnDum,
                     values_from = c(mean, sd)) %>% 
  dplyr::select(skim_variable, mean_1, sd_1, mean_0, sd_0) %>% 
  dplyr::mutate(skim_variable = c('救急', '機能', '支援', 
                                  '評価', '神経内科', '神経外科', 
                                  '病床数', '病床0', '大学病院')) %>% 
  kableExtra::kable(format = "html", 
                    booktabs = TRUE,
                    caption = "記述統計", 
                    col.names = c('変数', '平均', '標準偏差',
                                  '平均', '標準偏差'),
                    digits = 3) %>% 
  kableExtra::add_header_above(
    c(" ", "MRIを購買した病院" = 2, 
      "MRIを購買していない病院" = 2)) %>%
  kableExtra::kable_styling()

記述統計
	MRIを購買した病院		MRIを購買していない病院
変数	平均	標準偏差	平均	標準偏差
救急	0.712	0.453	0.302	0.459
機能	0.020	0.140	0.003	0.050
支援	0.048	0.213	0.004	0.063
評価	0.433	0.496	0.169	0.375
神経内科	0.393	0.489	0.122	0.327
神経外科	0.570	0.495	0.082	0.275
病床数	2.077	1.941	0.392	0.711
病床0	0.058	0.234	0.453	0.498
大学病院	0.042	0.200	0.005	0.067

各カテゴリーの病院のMRI保有割合

# 以下の表を作成するための関数
per_MRI_cal <- function(df){
  df %>% 
    dplyr::summarise(Total = n(),
                     MRIHos = sum(MRIOwnDum),
                     PerMRI = round(MRIHos / Total * 100, 2))
}

# 病院のカテゴリー毎に、総病院数、MRI所有病院数、その割合を確認
list(data_cleaned,
     data_cleaned %>% dplyr::filter(Kyukyu == 1),
     data_cleaned %>% dplyr::filter(Sien == 1),
     data_cleaned %>% dplyr::filter(Hyoka == 1),
     data_cleaned %>% dplyr::filter(DepNeurology == 1),
     data_cleaned %>% dplyr::filter(DepNeurosurgery == 1),
     # 元データの上位25%が1.2
     data_cleaned %>% dplyr::filter(LogNumBeds >= log(1.2)),
     data_cleaned %>% dplyr::filter(DaigakuDum == 1),
     data_cleaned %>% dplyr::filter(ZeroBedDum == 1)
     ) %>%
  purrr::map_dfr(.f =  per_MRI_cal) %>% 
  dplyr::mutate(Category = 
                  c('All', 'Kyukyu', 'Sien', 
                    'Hyoka', 'DepNeurology', 'DepNeurosurgery', 
                    'Top25%_NumBeds', 'DaigakuDum', 'ZeroBedDum')) %>% 
  dplyr::relocate(Category, .before = Total) %>% 
  kableExtra::kable(format = "html", 
                    booktabs = TRUE) %>% 
  kableExtra::kable_styling()

Category	Total	MRIHos	PerMRI
All	8862	3365	37.97
Kyukyu	4051	2392	59.05
Sien	182	160	87.91
Hyoka	2386	1456	61.02
DepNeurology	1987	1318	66.33
DepNeurosurgery	2365	1914	80.93
Top25%_NumBeds	2264	1932	85.34
DaigakuDum	165	140	84.85
ZeroBedDum	2674	196	7.33

5 Berry (1992)による推定

5.1 推定のためのデータ加工

Berry (1992) の推定を行うために、以下のような手順でデータを加工する。
- Step 0: 欠損のあるデータを削除
- Step 1: シミュレーションに必要な変数を作成、参入の順番に並び替え
- Step 2: 潜在的参入病院数を絞ってサブサンプル
- Step 3: シミュレーションを行いやすくするためにデータを縦に長く加工
Step 1, 3 については、それぞれの手順を関数として保存している。なぜなら、これらの手順は標準誤差を計算するためのブートストラップ法でも必要となるからである。
Step 2は本来不要だが、今回は推定にかかる時間を短縮するためにサンプルを絞っている。

# Step 0: 欠損のあるデータを削除
data_processed_pre <- 
  data_cleaned %>% 
  # 必要な変数のみ抽出
  dplyr::select(CityCode, Kyukyu, Kinou, Sien, Hyoka,
         DepNeurology, DepNeurosurgery, LogNumBeds, 
         ZeroBedDum, DaigakuDum,
         Menseki, LogPop, LogIncome, 
         MRIOwnDum) %>% 
  # 欠損を含む観察を削除
  na.omit()

# Step 1: シミュレーションに必要な変数を作成、参入の順番に並び替え
berry_process_data <- function(df){
  data_processed <- 
    df %>% 
    dplyr::arrange(CityCode) %>% 
    # 各病院に乱数を振る
    dplyr::mutate(TieEntryOrder = runif(nrow(df))) %>% 
    # 自治体毎の潜在的参入病院数、観察された参入病院数の変数を作成
    #（本誌のIm, Nmに対応）
    dplyr::group_by(CityCode) %>% 
    dplyr::mutate(NumPotenHos = as.numeric(n()),
                  NumEntryObs = sum(MRIOwnDum)) %>% 
    # データを参入の順番に並び替える
    # 病床数が多い順に並び替える。タイの場合は乱数を参照。
    dplyr::arrange(desc(LogNumBeds), desc(TieEntryOrder),
                 .by_group = TRUE) %>% 
    # incumbents move first の仮定の場合
    # dplyr::arrange(desc(MRIOwnDum), desc(LogNumBeds), desc(TieEntryOrder),
    #                .by_group = TRUE) %>% 
    # 自治体毎に、参入した順番に沿った病院idを振る
    dplyr::mutate(EntryOrderId = as.numeric(1:NumPotenHos)) %>% 
    dplyr::ungroup()
  return(data_processed)
}

# 以降の結果が変わらないように、乱数を固定
set.seed(1)
# Step 1で作成した関数の実行
data_processed <- berry_process_data(data_processed_pre)
# 不要なオブジェクトを削除
rm(data_processed_pre)

# Step 2: 潜在的参入病院数を絞ってサブサンプル
NumPotenHos_max <- 4
data_processed <- data_processed %>%
  # NumEntryObs_max 以下の観察に絞る
  dplyr::filter(EntryOrderId <= NumPotenHos_max) %>%
  # NumEntryObs, NumPotenHos の値を更新
  dplyr::group_by(CityCode) %>% 
  dplyr::mutate(NumEntryObs = sum(MRIOwnDum),
                NumPotenHos = as.numeric(n())) %>% 
  dplyr::ungroup()

各カテゴリーの病院のMRI保有割合（サブサンプル後のデータ）

# 病院のカテゴリー毎に、総病院数、MRI所有病院数、その割合を確認
list(data_processed,
     data_processed %>% dplyr::filter(Kyukyu == 1),
     data_processed %>% dplyr::filter(Sien == 1),
     data_processed %>% dplyr::filter(Hyoka == 1),
     data_processed %>% dplyr::filter(DepNeurology == 1),
     data_processed %>% dplyr::filter(DepNeurosurgery == 1),
     # 元データの上位25%が1.2
     data_processed %>% dplyr::filter(LogNumBeds >= log(1.2)),
     data_processed %>% dplyr::filter(DaigakuDum == 1),
     data_processed %>% dplyr::filter(ZeroBedDum == 1)
     ) %>%
  purrr::map_dfr(.f =  per_MRI_cal) %>% 
  dplyr::mutate(Category = 
                  c('All', 'Kyukyu', 'Sien', 
                    'Hyoka', 'DepNeurology', 'DepNeurosurgery', 
                    'Top25%_NumBeds', 'DaigakuDum', 'ZeroBedDum')) %>% 
  dplyr::relocate(Category, .before = Total) %>% 
  kableExtra::kable(format = "html", 
                    booktabs = TRUE) %>% 
  kableExtra::kable_styling()

Category	Total	MRIHos	PerMRI
All	4112	2248	54.67
Kyukyu	2477	1744	70.41
Sien	151	141	93.38
Hyoka	1426	1096	76.86
DepNeurology	1272	1015	79.80
DepNeurosurgery	1696	1469	86.62
Top25%_NumBeds	1857	1629	87.72
DaigakuDum	125	115	92.00
ZeroBedDum	646	43	6.66

# Step 3: シミュレーションを行いやすくするためにデータを縦に長く加工
# データを縦にns倍する
berry_expand_data <- function(df, ns){
  # 自治体の数
  M <- df %>% 
    dplyr::select(CityCode) %>% 
    dplyr::n_distinct() 
  # 病院の数（本誌のImのmについての和に対応）
  NumHos <- df %>%
    dplyr::distinct(CityCode, NumPotenHos) %>%
    dplyr::summarise(sum(NumPotenHos)) %>%
    as.numeric()
  
  # 誤差項をシミュレーション
  u_m0 <- rnorm(M * ns, mean = 0, sd = 1)
  u_mIm <- rnorm(NumHos * ns, mean = 0, sd = 1)
  u_df <- df %>% 
    # 自治体、病院idのみのdataframeを作成
    dplyr::select(CityCode, EntryOrderId) %>% 
    dplyr::group_by(CityCode) %>% 
    nest() %>% 
    # シミュレーション番号を振る（列数はns倍）
    tidyr::crossing(s = 1:ns) %>% 
    # 自治体毎の誤差項の変数を作成
    dplyr::mutate(u_m0) %>% 
    tidyr::unnest() %>% 
    # 病院毎の誤差項の変数を作成 
    dplyr::mutate(u_mIm)
  # 誤差項のdataframeを元のdataframeに統合
  data_expand <- df %>% 
    dplyr::left_join(u_df, by = c('CityCode', 'EntryOrderId'))
  return(data_expand)
}

# Step 3 で作成した関数の実行
ns <- 100
data_expand <- 
  berry_expand_data(
    data_processed, 
    ns)

5.2 Berry (1992) の目的関数の作成

# データとパラメターを所与としたBerry (1992) の目的関数の値を計算する関数を作成
berry_obj <- function(params, df){
  # パラメター
  alpha <- params[1:8]
  beta <- params[9:12]
  delta <- params[13]
  rho <- params[14]
  # 可変利潤、固定費用の作成（本誌のvとphiに対応）
  df <- df %>% 
    dplyr::mutate(
      var_profit = 
        beta[1] + beta[2] * Menseki + 
        beta[3] * LogPop + beta[4] * LogIncome +
        - delta * log(EntryOrderId) + rho * u_m0,
      fixed_cost = 
        alpha[1] * Kyukyu + 
        alpha[2] * Sien + alpha[3] * Hyoka +
        alpha[4] * DepNeurology + alpha[5] * DepNeurosurgery + 
        alpha[6] * LogNumBeds + alpha[7] * ZeroBedDum + 
        alpha[8] * DaigakuDum +
        - sqrt(1 - rho ^2) * u_mIm
      )
  # Berry (1992) による手法の目的関数を計算する
  obj <- df %>% 
    # 自治体、s、病院id、について、
    # 参入するかどうかの論理値型の列を作成
    dplyr::mutate(EntryDecision = var_profit > fixed_cost) %>% 
    # 自治体、s、について、不等式を満たす企業数を取得（本誌のN*と対応）
    # EntryDecision == FALSEという条件下で、EntryOrderIdの最小値はN*+1に対応するので、
    # 以下のような処理でN*+1を抽出している。論理値の条件分岐は処理速度が遅いため、
    # それを回避するために以下のようなコードになっている。
    # 潜在的参入企業数が10^5以上の場合は正しく動作しないので注意。
    dplyr::mutate(n_cand = (EntryDecision * 10^5 + 1) * EntryOrderId) %>%
    dplyr::group_by(CityCode, s) %>%
    dplyr::mutate(n_star = min(n_cand) - 1) %>%
    ungroup() %>%
    # 例外処理: 全て(EntryDecision == TRUE)となる場合は、
    # n_starをNumPotenHosに置き換える
    dplyr::mutate(n_star = ifelse(n_star >= 10^5, NumPotenHos, n_star)) %>% 
    # 自治体について、期待病院数を計算
    dplyr::group_by(CityCode) %>% 
    dplyr::summarise(n_exp = mean(n_star)) %>% 
    dplyr::left_join(df %>% dplyr::distinct(CityCode, NumEntryObs), 
              by = 'CityCode') %>% 
    # 実際の参入病院数と比較し、誤差の二乗を計算
    dplyr::mutate(dif = (NumEntryObs - n_exp)^2) %>% 
    dplyr::summarise(mean(dif)) %>% 
    as.numeric()
  return(obj)
}

5.3 推定

# 係数の名前
est_names <- 
  c('Kyukyu', 'Sien', 'Hyoka',
    'DepNeurology', 'DepNeurosurgery', 'LogNumBeds',
    'ZeroBedDum', 'DaigakuDum', 'Constant',
    'Menseki', 'LogPop', 'LogIncome',
    'delta', 'rho')

# 初期値
param_init <- 
  c(-0.612340533,   -5.525423772,   -0.505275676,   
    -0.32531026,    -1.04162392,    -0.991878025,   
    -3.87040966,    -1.272714254,   2.684741676,    
    0.040555764,    0.426448612,    -1.399627382,   
    0.990975782,    0.958075433)
# Obj: 0.372646328

以上の初期値の元で推定を行う。

tic.clearlog() 
tic()
berry_result <- 
  optim(params <- param_init,
        berry_obj, 
        method = "Nelder-Mead",
        df = data_expand)
# rhoの符号を正しくする
berry_result$par <- berry_result$par %>% 
  setNames(est_names)
berry_result$par[['rho']] <- abs(berry_result$par[['rho']])
berry_result$par <- berry_result$par %>% as.numeric()
toc()

## 243.81 sec elapsed

5.4 ブートストラップ法による標準誤差の計算

並列計算によってブートストラップを行う。

留意点： このブートストラップ計算には比較的時間がかかる点に留意されたい。並列計算に用いるコア数やCPUによっても異なるが、筆者のラップトップ(インテル Core i7-8650U、並列計算で６コア使用)の場合約３時間、ワークステーション(インテル Core i9-10980XE、並列計算で20コア使用)の場合約40分ほどかかる。計算時間や環境などに難がある方は、以下の2つのコードチャンクをスキップした上で、再現パッケージに含まれているBootStrap_Berry_Result.rdsファイルを用いて、以下のコードを読み進めてもらえれば幸いである。なお、RMDファイルにおける該当の2つのコードチャンクはeval = FALSEという設定にしているため、全体を走らせた場合には評価されない(コードとして実行されない)点に注意されたい。

最初のチャンクでは、並列計算のための設定を行う。ここで変数coresにおいて、並列計算に使うCPUコア数を設定する。なお、各ブートストラップサンプルにおける推定において乱数を利用するため、前と同様に乱数を固定する必要がある。しかしながら、並列計算の場合にはset.seed()では乱数を固定することができない。ここではdoRNGパッケージにおけるregisterDoRNG()を利用する。

# 並列計算のための設定を行う。
cores <- 20 #並列計算に使うCPUコアの個数。
registerDoParallel(cores)
registerDoRNG(12345)

並列計算によりブートストラップを行い、結果を保存する。並列計算にはforeachを使う。

# ブートストラップの回数
iter <- 100
# 自治体の数
M <- data_processed %>%
  dplyr::select(CityCode) %>%
  dplyr::n_distinct()
tic()
# 繰り返し推定を行う
bs_berry_result <-
  foreach::foreach (
    i = 1:iter,
    .export = ls(envir=parent.frame()),
      .packages = 'tidyverse'
    ) %dopar% {
      # リサンプルを行う
      bs_data <- data_processed %>%
        dplyr::distinct(CityCode) %>%
        dplyr::slice_sample(n = M, replace = TRUE) %>%
        dplyr::arrange(CityCode) %>%
        dplyr::mutate(id = 1:M) %>%
        dplyr::left_join(data_processed, by = 'CityCode') %>%
        # 推定ではCityCodeにより自治体を判別しているため名前を修正
        dplyr::rename(CityCodeOrigin = CityCode, CityCode = id)
      # リサンプルしたデータを推定のために加工
      bs_data_processed <-
        berry_process_data(bs_data)
      # 先程と同様に推定のためにデータを縦に広げる
      bs_data_expand <-
        berry_expand_data(
          bs_data_processed,
          ns)
      # 先程と同様に最適化により推定
      bs_berry_est <-
        optim(params <- param_init,
              berry_obj,
              method = "Nelder-Mead",
              df = bs_data_expand
              )$par
      # rhoの符号を正しくする
      bs_berry_est <- bs_berry_est %>%
        setNames(est_names)
      bs_berry_est[['rho']] <- abs(bs_berry_est[['rho']])
      return(bs_berry_est)
  }
toc()

# 結果を保存する。
saveRDS(bs_berry_result, file = "BootStrap_Berry_Result.rds")

ブートストラップの結果に基づいて標準誤差を計算する。

bs_berry_result <- readRDS("BootStrap_Berry_Result.rds")

# 標準誤差を計算
berry_se <- bs_berry_result %>%
  map_dfr(~ unlist(as.data.frame(.))) %>%
  summarise(across(everything(), sd)) %>%
  as.numeric()

6 推定結果

6.1 推定結果の確認

# 推定結果
tibble(
  name = c(est_names, 'Value of obj'),
  BR_est = c(berry_result$par, berry_result$value),
  # BR_se = berry_result_1129$BR_se
  BR_se = c(berry_se, NA)
  ) %>% 
  dplyr::mutate(name = c('救急', '支援', '評価', '神経内科', '神経外科', 
                         'log(病床数)', '病床0', '大学病院',
                         '定数項', '面積', 'log(人口)', 'log(課税所得)',
                         'delta', 'rho', '目的関数の値')) %>% 
  kableExtra::kable(format = "html", 
                    booktabs = TRUE,
                    caption = "Berry (1992) による推定結果", 
                    col.names = c('変数', '推定値', '標準誤差'),
                    digits = 3) %>%
  kable_styling()

Berry (1992) による推定結果
変数	推定値	標準誤差
救急	-0.589	0.065
支援	-5.588	0.125
評価	-0.497	0.083
神経内科	-0.349	0.124
神経外科	-1.048	0.105
log(病床数)	-0.969	0.042
病床0	-3.230	0.161
大学病院	-1.243	0.172
定数項	2.686	0.045
面積	0.033	0.014
log(人口)	0.418	0.017
log(課税所得)	-1.434	0.030
delta	0.965	0.060
rho	0.966	0.033
目的関数の値	0.370	NA

6.2 モデルによる予測の確認

# 各市場において病院の参入をシミュレートする関数
entry_sim <- function(df, delta){
  N <- df$NumPotenHos %>% unique()
  prof_vec <- df$prof_excl_comp
  EntryCond <- c()
  n <- 1
  for (i in 1:N) {
    net_prof <- prof_vec[i] - delta * log(n)
    if (net_prof > 0) { 
      # 参入する場合
      EntryCond <- c(EntryCond, 1)
      n <- n + 1
    } else { 
      # 参入しない場合
      EntryCond <- c(EntryCond, 0)
    }
  }
  return(EntryCond)
}

# 推定したパラメターを取得
berry_est <- berry_result$par
alpha <- berry_est[1:8]
beta <- berry_est[9:12]
delta <- berry_est[13]
rho <- berry_est[14]

# モデルによる参入予測確率を含めたdataframeを作成
data_predicted <- data_expand %>% 
  # 可変利潤、固定費用の作成（本誌のvとphiに対応）
  dplyr::mutate(
    var_profit = 
      beta[1] + beta[2] * Menseki + 
      beta[3] * LogPop + beta[4] * LogIncome +
      rho * u_m0,
    fixed_cost = 
      alpha[1] * Kyukyu  + 
      alpha[2] * Sien + alpha[3] * Hyoka +
      alpha[4] * DepNeurology + alpha[5] * DepNeurosurgery + 
      alpha[6] * LogNumBeds + alpha[7] * ZeroBedDum + 
      alpha[8] * DaigakuDum +
      - sqrt(1 - rho ^2) * u_mIm,
    # 競争の効果を除いた利潤を計算
    prof_excl_comp = var_profit - fixed_cost
    ) %>% 
  # entry_sim を用いて、自治体、s毎に、企業の参入行動をシミュレート
  dplyr::group_by(CityCode, s) %>% 
  tidyr::nest() %>% 
  dplyr::mutate(EntryCond = purrr::map(data, entry_sim, delta = delta)) %>% 
  tidyr::unnest(cols = c(data, EntryCond)) %>% 
  dplyr::ungroup() %>% 
  # 参入の予測確率を計算
  dplyr::group_by(CityCode, EntryOrderId) %>% 
  dplyr::mutate(EntryProb = mean(EntryCond)) %>% 
  dplyr::ungroup() %>% 
  dplyr::mutate(EntryPred = ifelse(EntryProb >= 0.5 ,1 ,0)) %>% 
  dplyr::select(-c(s, u_m0, u_mIm, var_profit, fixed_cost, prof_excl_comp, EntryCond)) %>% 
  dplyr::distinct(across(everything()))

# 参入病院数毎に、自治体がいくつ観察されるかを集計
data_predicted %>% 
  dplyr::group_by(CityCode) %>% 
  dplyr::summarise(Actual = sum(MRIOwnDum),
                   Predict = sum(EntryPred)) %>% 
  tidyr::pivot_longer(cols = c('Actual', 'Predict')) %>% 
  ggplot2::ggplot(aes(x = value, fill = name)) + 
  ggplot2::geom_histogram(binwidth = 0.5,
                          position = "dodge") + 
  coord_flip() + 
  scale_x_reverse()

# 以下の表を作成するための関数
actual_pred_cal <- function(df){
  df %>% 
    dplyr::summarise(Actual = sum(MRIOwnDum),
                     Predict = sum(EntryPred))
}

# 病院のタイプ別の観察された参入病院数と予測された参入病院数
list(data_predicted,
     data_predicted %>% dplyr::filter(Kyukyu == 1),
     data_predicted %>% dplyr::filter(Sien == 1),
     data_predicted %>% dplyr::filter(Hyoka == 1),
     data_predicted %>% dplyr::filter(DepNeurology == 1),
     data_predicted %>% dplyr::filter(DepNeurosurgery == 1),
     # 元データの上位25%が1.2
     data_predicted %>% dplyr::filter(LogNumBeds >= log(1.2)),
     data_predicted %>% dplyr::filter(DaigakuDum == 1),
     data_predicted %>% dplyr::filter(ZeroBedDum == 1)
     ) %>%
  purrr::map_dfr(.f =  actual_pred_cal) %>% 
  dplyr::mutate(Category = 
                  c('All', 'Kyukyu', 'Sien', 
                    'Hyoka', 'DepNeurology', 'DepNeurosurgery', 
                    'Top25%_NumBeds', 'DaigakuDum', 'ZeroBedDum')) %>% 
  tidyr::pivot_longer(cols = c('Actual', 'Predict')) %>%
  ggplot2::ggplot(aes(x = Category, y = value, fill = name)) + 
  geom_bar(position="dodge", stat = "identity") + 
  coord_flip()

  # theme_minimal()

# 特定の病院について平均二乗誤差を表示
tibble(
  category = c('Total', 'DepNeurology', 'DepNeurosurgery', 'ZeroBedDum', 'DaigakuDum'),
  mse = c(data_predicted %>% 
            dplyr::summarise(mean((MRIOwnDum - EntryProb) ^ 2)) %>% 
            as.numeric() %>% round(3),
          data_predicted %>% 
            dplyr::filter(DepNeurology == 1) %>% 
            dplyr::summarise(mean((MRIOwnDum - EntryProb) ^ 2)) %>% 
            as.numeric() %>% round(3),
          data_predicted %>% 
            dplyr::filter(DepNeurosurgery == 1) %>% 
            dplyr::summarise(mean((MRIOwnDum - EntryProb) ^ 2)) %>% 
            as.numeric() %>% round(3),
          data_predicted %>% 
            dplyr::filter(ZeroBedDum == 1) %>% 
            dplyr::summarise(mean((MRIOwnDum - EntryProb) ^ 2)) %>% 
            as.numeric() %>% round(3),
          data_predicted %>% 
            dplyr::filter(DaigakuDum == 1) %>% 
            dplyr::summarise(mean((MRIOwnDum - EntryProb) ^ 2)) %>% 
            as.numeric() %>% round(3)
          )
)

7 反実仮想分析

# 神経内科と神経外科を持つ病院はMRIの所有を義務付けられたとする
# MRI所有義務化病院が既に参入しているとして、他の病院が参入するかどうかをシミュレート
data_cf <- data_expand %>% 
  dplyr::group_by(CityCode, s) %>% 
  # MRI所有義務化病院が先になるように並び替える
  dplyr::arrange(desc(DepNeurology), desc(DepNeurosurgery), 
                 desc(LogNumBeds), desc(TieEntryOrder),
                 .by_group = TRUE) %>% 
  # 参入する順番を更新する
  dplyr::mutate(EntryOrderId = as.numeric(1:NumPotenHos)) %>% 
  dplyr::ungroup() %>% 
  # 可変利潤、固定費用の作成（本誌のvとphiに対応）
  dplyr::mutate(
    var_profit = 
      beta[1] + beta[2] * Menseki + 
      beta[3] * LogPop + beta[4] * LogIncome +
      rho * u_m0,
    fixed_cost = 
      alpha[1] * Kyukyu +  
      alpha[2] * Sien + alpha[3] * Hyoka +
      alpha[4] * DepNeurology + alpha[5] * DepNeurosurgery + 
      alpha[6] * LogNumBeds + alpha[7] * ZeroBedDum + 
      alpha[8] * DaigakuDum +
      - sqrt(1 - rho ^2) * u_mIm,
    # 競争の効果を除いた利潤を計算
    prof_excl_comp = var_profit - fixed_cost
    ) %>% 
  # MRI所有義務化病院が必ず参入するために大学病院の利潤を大きな値に置き換え
  dplyr::mutate(prof_excl_comp = 
                  ifelse(DepNeurology == 1 | DepNeurosurgery == 1
                         , 10^5, prof_excl_comp)) %>% 
  # entry_sim を用いて、自治体、s毎に、企業の参入行動をシミュレート
  dplyr::group_by(CityCode, s) %>% 
  tidyr::nest() %>% 
  dplyr::mutate(EntryCond = purrr::map(data, entry_sim, delta = delta)) %>% 
  tidyr::unnest(cols = c(data, EntryCond)) %>% 
  dplyr::ungroup() %>% 
  # 参入の予測確率を計算
  dplyr::group_by(CityCode, EntryOrderId) %>% 
  dplyr::mutate(EntryProb = mean(EntryCond),
                EntryPred = ifelse(EntryProb >= 0.5 ,1 ,0)) %>% 
  dplyr::ungroup() %>% 
  dplyr::select(-c(s, u_m0, u_mIm, var_profit, fixed_cost, prof_excl_comp, EntryCond)) %>% 
  dplyr::distinct(across(everything()))

# 以下の表を作成するための関数
per_decline_cal <- function(df){
  df %>% 
    dplyr::summarise(
      Total = n(),
      Actual = sum(MRIOwnDum),
      Fit_Pred = sum(Fit_Pred),
      CF_Pred = sum(EntryPred),
      CF_Pred_am_actual = sum(MRIOwnDum * EntryPred)
      )
}

data_cf <- 
  data_cf %>% 
  # 元のモデルでの参入の予測を表す変数を加える
  dplyr::left_join(data_predicted %>% 
                     dplyr::select(CityCode, TieEntryOrder, EntryPred) %>% 
                     dplyr::rename(Fit_Pred = EntryPred),
                   by = c("CityCode", "TieEntryOrder")) %>%
  # 政策の対象かどうかを表すダミー
  dplyr::mutate(PolicyTarget = 
                  ifelse(
                    (DepNeurology == 1 | DepNeurosurgery == 1) & Fit_Pred != 1,
                      #DepNeurology == 1 | DepNeurosurgery == 1) & MRIOwnDum != 1,
                    1, 0)
                )

# 病院のカテゴリー毎に、総病院数(Total)、
# 元のモデルにより予測されたMRI所有病院数(Fit_Pred)、
# 反実仮想分析により予測されたMRI所有病院数(CF_Pred)、
# MRI所有病院数の変化率(PerChange)、を確認
data_cf_sub <- data_cf %>% 
  dplyr::filter(PolicyTarget == 0)

list(data_cf_sub,
     data_cf_sub %>% dplyr::filter(Kyukyu == 1),
     data_cf_sub %>% dplyr::filter(Sien == 1),
     data_cf_sub %>% dplyr::filter(Hyoka == 1),
     data_cf_sub %>% dplyr::filter(DepNeurology == 1),
     data_cf_sub %>% dplyr::filter(DepNeurosurgery == 1),
     data_cf_sub %>% dplyr::filter(!(DepNeurology == 1 | DepNeurosurgery == 1)),
     # 元データの上位25%が1.2
     data_cf_sub %>% dplyr::filter(LogNumBeds >= log(1.2)),
     data_cf_sub %>% dplyr::filter(DaigakuDum == 1),
     data_cf_sub %>% dplyr::filter(ZeroBedDum == 1)
     ) %>%
  purrr::map_dfr(.f =  per_decline_cal) %>% 
  dplyr::mutate(Category = 
                  c('All', 'Kyukyu', 'Sien', 'Hyoka',
                    'DepNeurology', 'DepNeurosurgery', 
                    'Not_(DepNeurology_and_DepNeurosurgery)',
                    'Top25%_NumBeds', 'DaigakuDum', 'ZeroBedDum')) %>% 
  dplyr::select(Category, Fit_Pred, CF_Pred) %>% 
  dplyr::mutate(PerChange = round((CF_Pred/Fit_Pred - 1) * 100, 2)) %>%
  kableExtra::kable(format = "html",
                    booktabs = TRUE,
                    digits = 3) %>%
  kableExtra::kable_styling()

Category	Fit_Pred	CF_Pred	PerChange
All	2301	2255	-2.00
Kyukyu	1887	1861	-1.38
Sien	151	151	0.00
Hyoka	1201	1187	-1.17
DepNeurology	1076	1076	0.00
DepNeurosurgery	1619	1619	0.00
Not_(DepNeurology_and_DepNeurosurgery)	429	383	-10.72
Top25%_NumBeds	1790	1769	-1.17
DaigakuDum	123	123	0.00
ZeroBedDum	3	3	0.00

経済セミナー連載「実証ビジネス・エコノミクス」
第６回「戦略は企業の特性が決める：参入ゲームの推定［応用編］」

上武康亮・遠山祐太・若森直樹・渡辺安虎

最終更新: 2022-01-19

1 はじめに

1.1 留意事項

1.2 全体の流れ

1.3 謝辞

2 Rに関する下準備

3 データの準備

4 企業レベル変数の記述統計

5 Berry (1992)による推定

5.1 推定のためのデータ加工

5.2 Berry (1992) の目的関数の作成

5.3 推定

5.4 ブートストラップ法による標準誤差の計算

6 推定結果

6.1 推定結果の確認

6.2 モデルによる予測の確認

7 反実仮想分析

経済セミナー 連載「実証ビジネス・エコノミクス」 第６回 「戦略は企業の特性が決める：参入ゲームの推定［応用編］」

上武康亮・遠山祐太・若森直樹・渡辺安虎

最終更新: 2022-01-19

1 はじめに

1.1 留意事項

1.2 全体の流れ

1.3 謝辞

2 Rに関する下準備

3 データの準備

4 企業レベル変数の記述統計

5 Berry (1992)による推定

5.1 推定のためのデータ加工

5.2 Berry (1992) の目的関数の作成

5.3 推定

5.4 ブートストラップ法による標準誤差の計算

6 推定結果

6.1 推定結果の確認

6.2 モデルによる予測の確認

7 反実仮想分析

経済セミナー連載「実証ビジネス・エコノミクス」
第６回「戦略は企業の特性が決める：参入ゲームの推定［応用編］」