class: title-slide-section, center, middle
name: logistics
# イントロダクション


---
## イントロダクション

- これまで：統計的因果推論に基づく分析
    - 回帰分析、自然実験を用いて因果効果・政策効果を推定する。
    

- これから：構造推定アプローチ
    - ミクロ経済学モデルのパラメタの推定及び反実仮想シミュレーション
    

- 最終的なゴール：両アプローチの強み・弱み、そして使い分けを学ぶ。

---
## 非線形モデルの推定と最尤法

- ミクロ経済理論から得られるモデルはしばしば非線形になる。
    - ある被説明変数を説明変数の線形関数で書けるとは限らない。
    
    
- そのような場合に、しばしば最尤法によるパラメタ推定を行う。


- 授業の流れ：
    1. 最尤法と簡単な例
    2. 二項選択モデル：ロジットモデルとプロビットモデル
    3. Rプログラミング：政治的選好
    4. (時間と応相談)ケース：テック企業の価格実験


---
class: title-slide-section, center, middle
name: logistics
# 最尤法の入門

---
## セットアップ

-   確率変数 $y$ を考える。その密度関数を $f(y;\theta)$ とし、 $\theta$ は未知パラメタ。


-   密度関数は正しく定式化されていると仮定する。


-   データ $\left\{ y_{i}\right\} _{i=1}^{n}$ が得られ、各観測は独立かつ同一 (i.i.d.) に $f(y;\theta_{0})$ から生成されている。
    -   $\theta_{0}$ は未知のパラメタ


-   ゴール: 未知パラメタ $\theta_{0}$ を データ $\left\{ y_{i}\right\} _{i=1}^{n}$ を用いて推定する。

---
## 最尤推定量 (Maximum Likelihood Estimator, MLE)

-   データ $\left\{ y_{i}\right\} _{i=1}^{n}$ に関する尤度関数を以下のように定義する。
$$\begin{aligned}
    L(\theta) & =\prod_{i=1}^{n}f(y_{i};\theta),\end{aligned}$$

    -   データ
        $\left\{ y_{i}\right\} _{i=1}^{n}$ を観測する確率を パラメタ $\theta$ の関数として表したもの。

-   対数を取った**対数尤度関数**を用いることが多い。
    $$\log L(\theta)=\sum_{i=1}^{n}\log f(y_{i};\theta).$$

-   最尤推定量を以下のように定義する。
    $$\hat{\theta}_{n}=\operatorname{argmax}_{\theta\in\Theta}\log L(\theta).$$

---
## 例:コイン投げ

-   以下の確率変数 $Y_{i}$ を考える。 
$$Y_{i}=\begin{cases}
    1 & with\ prob\ p\\
    0 & with\ prob\ 1-p
    \end{cases}$$

-   データとして $N$ 回のコイン投げの結果がわかるとしよう。 $\{y_{i}\}_{i=1}^{N}$

-   尤度関数は 
$$L(\theta) = \prod_{i=1}^{n}P(Y=y_{i};p)=\prod_{i=1}^{n}p^{\mathrm{1}\{y_{i}=1\}}(1-p)^{\mathrm{1}\{y_{i}=0\}}$$ 

---
## コイン投げモデルにおける対数尤度

- 対数尤度関数は
$$\begin{eqnarray} 
\log L(\theta) &=& \sum_{i=1}^{n}\log\left\{ p^{\mathrm{1}\{y_{i}=1\}}(1-p)^{\mathrm{1}\{y_{i}=0\}}\right\} \\ &=& \sum_{i=1}^{n}\mathrm{1}\{y_{i}=1\}\log p+\mathrm{1}\{y_{i}=0\}\log(1-p) \\ 
&=& N_{1}\log p+N_{0}\log(1-p)
\end{eqnarray}$$
ここで、 $N_{1}=\sum_{i=1}^{N}\mathbf{1}\{y_{i}=1\},N_{0}=\sum_{i=1}^{N}\mathbf{1}\{y_{i}=0\}$

- 対数尤度を最大にする $\hat{p}_{MLE}$ は $N_1 / N$ となる。



---
## 最尤推定量の性質

-   最尤推定量の定義：
    $$\hat{\theta}=\operatorname{argmax}_{\theta\in\Theta}\log L(\theta|x_{i}).$$


-   最尤推定量も一定の仮定のもとで、以下のような性質を満たす。
    -   一致性 (Consistency): As the sample size $N$ gets larger, the estimator
        is getting closer to the true value $\theta_{0}$
    -   漸近正規性 (Asymptotic normality): The distribution of $\hat{\theta}$ is
        approximated by a normal distribution as the sample size gets
        larger.


- 詳しくは補足資料及び各種教科書を参照。

---
class: title-slide-section, center, middle
name: logistics
# 二項被説明変数のモデル


---
## 制限的従属変数 (limited dependent variable)

- 最尤推定量がよく用いられるのは、従属変数が制限的なケース


- 例１：二値(ダミー)変数、カテゴリカル変数
    - 二項選択モデル、多項選択モデル、順序付き選択モデル


- 例２：打ち切り変数
    - 働いていない場合、賃金はゼロと記録される。
    - トービットモデル、ヘックマンのサンプルセレクションモデル


- 例３：整数値のみを取る値 
    - ポアソン回帰モデル
    
    
- 本コースでは、二項選択・多項選択モデルを中心に扱う。
    - その他については、Wooldridge "Introductory Econometrics"を参照。

---
## 被説明変数が二値変数(ダミー変数)のモデル

-   被説明変数がダミー変数 (0 or 1) となるモデルを考えていく。


- 例１：働くか働かないか


- 例２：ある製品を購入するかしないか。


- 例３：処置（トリートメント）の有無。
  - 傾向スコア $P(D=1 | X=x)$ の推定においては二項選択モデルがしばしば用いられる。
  
  
- 例４：ローン審査。金融機関がローンを承認するか否か。



---
## 線形確率モデル

-   被説明変数がダミー変数でも線形回帰モデルを使うことも可能。  $$y_{i}=\beta_0 + \beta_1 x_{1i} + \cdots + \beta_K x_{Ki} +\epsilon_{i}$$
    -   **線形確率モデル**と呼ばれる。


-   問題点：モデルが予測する確率が $[0,1]$ 区間の外に出る可能性。
  - 購買行動の予測や傾向スコアの計算には不適


-   とはいえ、線形確率モデルはしばしば用いられる。
    - 利点：操作変数を使った内生性の処理や、ダミー変数（固定効果）を柔軟に扱える。
    - もし特定の変数の係数(因果効果)のみに興味があるならば、これで十分なことも多い。


---
## 二値従属変数のモデル

-   以下のモデルを考えよう
    $$d_{i}=\boldsymbol{1}\{\beta_0 + \beta_1 x_{1i} + \cdots + \beta_K x_{Ki}+\epsilon_{i}\geq0\}.$$ 
    ここで、 $\boldsymbol{1}\{\cdot\}$ は指示関数 (1 if TRUE and 0
    otherwise)

- 解釈1: 被説明変数 $d_{i}$ が上記のルールを満たすときに1となる。


- 解釈2: $d=1$ という「選択肢」を取ったときの効用が $x_{i}'\beta+\epsilon_{i}$, $d=0$のときには効用が $0$ となるもとで、自分の効用が最大化するような選択を行っている。
    - 来週、離散選択モデルにおいて詳しく。



---
## 誤差項に関する仮定

- このモデルはパラメタに関して非線形になっている。OLSでは推定不可。


- 誤差項 $\epsilon$ に分布の仮定を置くことで、最尤法で推定を行う。


- 誤差項  $\epsilon$ の分布関数を $G(x) = P(\epsilon \leq x)$ と定義


- また、 $x_{i}'\beta = \beta_0 + \beta_1 x_{1i} + \cdots + \beta_K x_{Ki}$ と定義。


- $d=1$ となる確率は

$$
\begin{eqnarray}
P(d = 1 | x) &=& P(x_{i}'\beta+\epsilon_{i} > 0 |x) \\ 
&=& P( \epsilon_{i} > - x_{i}'\beta | x) \\
&=& 1-G(-x_{i}'\beta) \\
&=& G(x_{i}'\beta)
\end{eqnarray}
$$
    
---
## 代表的なモデル


-   プロビットモデル:  $\epsilon$ が標準正規分布に従う。分布関数は $$G(x) = \Phi (x) = \int_{- \infty}^x \phi(z) dz$$ ここで $\phi(x)$ は標準正規分布の密度関数

-   ロジットモデル: $\epsilon$ がロジスティック分布に従う。分布関数は $$G(x)=P(\epsilon\leq x)=\frac{\exp(z)}{1+\exp(z)}$$


---
## プロビットモデル

-   仮定： $\epsilon\sim N(0,1)$.

-   このとき、 $d_{i}=1$ となる確率は
$$\begin{aligned}
    P(d_{i}=1) & =P(\epsilon\geq-x_{i}'\beta)\\
     & =1-\Phi\left( -x_{i}'\beta\right)\\
     & =\Phi\left(x_{i}'\beta \right).\end{aligned}$$
    ここで $\Phi(\cdot)$ は標準正規分布の分布関数


---
## ロジットモデル

-   仮定 $\epsilon_{i}$ がロジスティック分布に従う。

-   このとき、 $d_{i}=1$ となる確率は
    $$P(d_{i}=1)=\frac{\exp\left( x_{i}'\beta \right)}{1+\exp\left( x_{i}'\beta \right)}.$$

---
## パラメタの解釈・限界効果

- 線形モデルのケースと異なって、パラメタの値そのものは左辺への直接的な影響ではない。


- 確率を $P(d_i = 1 | x_i ) = G(x_{i}'\beta)$ と書く。


- 変数 $x_{ki}$ の確率への効果は
$$\frac{\partial P(d_i = 1 | x_i )}{\partial x_{ki}} = G'(x_{i}'\beta)\times \beta_k$$
  - プロビット： $G'(x_{i}'\beta) = \phi(x_{i}'\beta)$ (pdf of normal)
  - ロジット： $G'(x_{i}'\beta) = \frac{\exp\left( x_{i}'\beta \right)}{1+\exp\left( x_{i}'\beta \right)} \left( 1- \frac{\exp\left( x_{i}'\beta \right)}{1+\exp\left( x_{i}'\beta \right)} \right)$


---
## 限界効果の計算

- 個々の観測ごとに計算し、それの平均をとる
$$\frac{1}{N} \sum_{i=1}^N \frac{\partial P(d_i = 1 | x_i )}{\partial x_{ki}} = \frac{1}{N} \sum_{i=1}^N G'(x_{i}'\beta)\times \beta_k$$

- サンプルの平均 $\bar{x}$ で評価する
$$\frac{\partial P(d_i = 1 | x_i = \bar{x} )}{\partial x_{ki}} =  G'(\bar{x}'\beta)\times \beta_k$$



---
class: title-slide-section, center, middle
name: logistics
# 二項被説明変数のモデルの推定


---
## 推定に先立ち：パラメタの正規化(ノーマリゼーション)


- モデルの性質上、識別できる・できないパラメタを区別する必要がある。


- プロビットにおいて $\epsilon\sim N(\mu,\sigma^{2})$ と考え、追加的なパラメタ $(\mu,\sigma)$ を考える。

- このとき、 $d_{i}=1$ となる確率は
$$\begin{aligned}
    P(d_{i}=1) & =P(\epsilon\geq-x_{i}'\beta)\\
     & =P\left(\frac{\epsilon - \mu}{\sigma}\geq\frac{-x_{i}'\beta - \mu}{\sigma}\right)\\
     & =\Phi\left( \frac{\beta_0 + \beta_1 x_{1i} + \cdots + \beta_K x_{Ki} + \mu }{\sigma}\right).\end{aligned}$$

---
## location and scale normalization

- 確率再掲

$$P(d_{i}=1) =\Phi\left( \frac{\beta_0 + \beta_1 x_{1i} + \cdots + \beta_K x_{Ki} + \mu }{\sigma}\right)$$

- モデル上、定数項 $\beta_0$ と 誤差項の平均 $\mu$ は別個に推定できない。
  - $\mu = 0$ という正規化が必要 (**location normalization**)
  
  
- モデル上、パラメタと誤差項の標準偏差の比率 $\beta / \sigma$ という形で現れる。
  - すべての $\beta$ と $\sigma$ を定数倍しても全くおなじになる。
  - $\sigma = 1$ という正規化が必要 (**scale normalization**)

---
## 二項選択モデルの尤度関数

-   データ： $\{d_{i},x_{i}\}_{i=1}^{N}$.

-   尤度関数
    $$L(\theta)=\prod_{i=1}^{N}P(d_{i}=0)^{\mathbf{1}\{d_{i}=0\}}P(d_{i}=1)^{\mathbf{1}\{d_{i}=1\}}$$

-   ロジットモデルの場合：
    $$L(\theta)=\prod_{i=1}^{N}\left(\frac{1}{1+\exp\left(x_{i}'\beta\right)}\right)^{\mathbf{1}\{d_{i}=0\}}\left(\frac{\exp\left(x_{i}'\beta\right)}{1+\exp\left(x_{i}'\beta\right)}\right)^{\mathbf{1}\{d_{i}=1\}}$$

---
## 対数尤度関数と数値最適化

- ロジットモデルの対数尤度関数

$$\log L(\theta) = \sum_{i=1}^N \mathbf{1}\{d_{i}=0\} \log \left(\frac{1}{1+\exp\left(x_{i}'\beta\right)}\right)+ \mathbf{1}\{d_{i}=1\}\left(\frac{\exp\left(x_{i}'\beta\right)}{1+\exp\left(x_{i}'\beta\right)}\right)$$

- この対数尤度関数を最大にするパラメタは解析的には得られない。-> 数値計算


- ロジットモデルでは対数尤度関数はパラメタに関して凹関数であり、数値計算は比較的容易。


---
## (おまけ) 二項選択モデルにおける内生性と操作変数


- 尤度関数を記述する上で、説明変数は外生(誤差と相関なし)を仮定している。


- 二項選択モデルにおける説明変数に内生性がある場合どうするか？


- コントロール関数アプローチ(control function)と呼ばれるものがしばしば使われる。
    - Step 1: 内生変数を操作変数と外生変数に回帰する。
    - Step 2: Step1で得られた"残差"を、元のモデルに追加的な説明変数として入れ、最尤法で推定する。
    
    
- 詳しくは 大学院レベルWooldridge などを参照。



---
class: title-slide-section, center, middle
name: logistics
# コード例：自民党支持の決定要因

---
## 例：自民党支持の決定要因

- 松浦「Rによるデータ分析入門」の第四章練習問題１より


- データファイル`political_party,csv` は以下よりダウンロード可
  - http://www.tokyo-tosho.co.jp/books/978-4-489-02424-5/
  - 履修者はMoodleを参照


- データ：東大社研若年者パネル調査の公開データ(恐らく疑似データ)
  - 申請すれば正式なデータが取得可能。


- 内容
  - 線形モデルの推定
  - ロジット・プロビットモデルの推定、係数比較
  - 限界効果の計算
  
  

---
## 変数一覧

| 変数名          | 説明       |
|--------------|----------|
| female       | 女性ダミー    |
| age          | 年齢       |
| income       | 所得（単位万円） |
| educ         | 教育年数     |
| value_family | 伝統的家族感   |
| social_class | 社会階層     |
| LDP          | 自民党支持ダミー |
| DPJ          | 民主党支持ダミー |
| COM          | 共産党支持ダミー |
| marriage     | 結婚ダミー    |


---
## 変数一覧の補足

- 伝統的家族観：「男性は収入を得て、女性は家庭と家族の面倒をみるべき」に対して「そう思う」が５、「そう思わない」が１を取る。

- 社会階層：「仮に社会全体を10段階で分けると、あなたはどこに位置すると思いますか？」に対して、10が一番上、1が一番下。


---
## コードの準備

- `fixest`パッケージにおける`feglm`関数を利用する。
  - Rに元から入っている`glm`関数のfixestバージョン
  - `glm`は generalized linear model (一般化線形モデル) の略称
- 限界効果の計算に `marginaleffects`パッケージを利用する。

```{r, eval = TRUE}
library(tidyverse)
library(fixest)
library(summarytools)
library(marginaleffects) 
dataf <- readr::read_csv("data/political_party.csv")
```


---
## 記述統計

```{r}
summarytools::descr(dataf, transpose = TRUE, stats = c("Mean", "sd", "Min", "Max"))
```

---
## 線形確率モデルで推定

- まずは線形確率モデル

```{r}
result1 <- feols(LDP ~ female + educ + marriage + income + 
                   age + value_family + social_class,
                 data=dataf, vcov = "hetero")
```

---
## 推定結果

```{r}
etable(result1)
```


---
## 線形確率モデルの予測確率

```{r}
pred_linprob = predict(result1)
summary(pred_linprob)
```


---
## ロジットモデル・プロビットモデル

```{r}
result2 <- feglm(　LDP ~ female + educ + marriage + income + 
                   age + value_family + social_class,
                  family=binomial(link=logit),
                 data=dataf, vcov = "hetero")

result3 <- feglm(　LDP ~ female + educ + marriage + income + 
                   age + value_family + social_class,
                  family=binomial(link=probit),
                 data=dataf, vcov = "hetero")

```

---
## 推定結果の比較

```{r}
etable(list(result1, result2, result3))
```

---
## ロジットモデル・プロビットモデルの予測確率

```{r}
pred_logit = predict(result2)
pred_probit = predict(result3)

summary(tibble( pred_linprob, pred_logit, pred_probit) )
```

---
## ロジットモデルの限界効果

- パッケージ `marginaleffects`を使おうとしたらエラーが出てうまくいかなかったので、手計算でやります。


- 変数 $x_{ki}$ の確率への効果は
$$
\begin{eqnarray}
\frac{\partial P(d_i = 1 | x_i )}{\partial x_{ki}} &=& G'(x_{i}'\beta)\times \beta_k \\
&=& \frac{\exp\left( x_{i}'\beta \right)}{1+\exp\left( x_{i}'\beta \right)} \left( 1- \frac{\exp\left( x_{i}'\beta \right)}{1+\exp\left( x_{i}'\beta \right)} \right) \times \beta_k \\
&=& P(d_i = 1 |x_i) \times \left( 1 - P(d_i = 1 |x_i) \right) \times \beta_k
\end{eqnarray}
$$
- よって、各観測における予測確率をまず計算する。

---
## 続・限界効果

- 各観測の限界効果の平均値を計算しよう。
$$\frac{1}{N} \sum_{i=1}^N \frac{\partial P(d_i = 1 | x_i )}{\partial x_{ki}} = \left( \frac{1}{N} \sum_{i=1}^N G'(x_{i}'\beta) \right) \times \beta_k$$


```{r}
dataf %>%
  mutate( pred_logit_prob = predict(result2)) -> dataf
avg_prob_1_prob = mean( dataf$pred_logit_prob * ( 1-dataf$pred_logit_prob ))

marg_effect = avg_prob_1_prob*result2$coefficients

```

---
## 比較：ロジットの限界効果と線形確率モデル

```{r}
comp <- cbind(marg_effect, result1$coefficients)
colnames(comp) <- c("logit", "linear prob")
print(comp)
```



---
class: title-slide-section, center, middle
name: logistics
# ケース：テック企業の価格実験

---
## ケース：オンライン求人サイトのプライシング

-   論文：Dube and Misra (2021) Personalized Pricing and Consumer Welfare


-   シカゴ大学の研究者たちがZiprecruiter.comというオンライン求人サイトと協力してプライシングの実験を行った。


-   流れ
    -   Ziprecruiter.com の紹介
    -   実験のデザイン
    -   実験の結果


---
## Ziprecruiter.com

-   オンライン求人サイトhttps://www.ziprecruiter.com/
    -   求職者(jobseeker): 無料
    -   求人企業(employers): 有料。分析の主眼


-   仕組み：月額サブスクリプション
    -   実験前：99ドル／月
    -   キャンセル自由
    -   約4万企業の有料会員
    -   限界費用はほぼゼロ(収入最大化)


---
## ページトップ

.center[
<img src = "figuretable/fig_zip_top.png" width="800">
]


---
## 求職者側のWebsite (２２年２月現在)

.center[
<img src = "figuretable/fig_zip_search.png" width="800">
]

---
## 2015年の実験

-   2015年春：プライシングの再検討を始める。
    -   2015年9月：実験１（需要推定）
    -   2015年11月：実験２（プライシングのValidation）
    

-   実験対象：
    -   Ziprecruiter.comに初めて加入
    -   Starter企業： 従業員50人未満。平均して1-3のポジションを募集


-   サブスクリプションの仕組み：
    -   ホームページ上に価格は出ていない。
    -   種々の情報の入力（仕事のタイプなど）する。
    -   支払画面において**初めて価格が提示**され、加入するならばクレジットカード情報を入力


---
## 実験１：ランダム価格

-   対照群・コントロールグループ： 99ドル
-   処置群・トリートメントグループ：19ドルから399ドルまでの価格でRandomizeする。

.center[
<img src = "figuretable/tab_experiment_price.png" width="350">
]



---
## 実験１の生結果：需要（左）と収入（右）


.center[
<img src = "figuretable/fig_raw_result.png" width="1000">
]



<!-- --- -->
<!-- ## (参考)2015年9月から12月にかけてのRetention (継続率) -->

<!-- .center[ -->
<!-- <img src = "figuretable/fig_acquisition.png" width="1000"> -->
<!-- ] -->


<!-- --- -->
<!-- ## (参考)2015年9月から12月の割引現在価値(継続率を考慮) -->


<!-- .center[ -->
<!-- <img src = "figuretable/fig_revenue.png" width="800"> -->
<!-- ] -->

---
## ランダム化比較試験から需要推定へ

- ここまで：ランダム化比較試験により、価格が加入率(需要)に与える影響がわかった。


- では、具体的にいくらの値段に設定するべきか？ -> モデルが必要！！


- 二項ロジットモデルを用いて需要関数を推定しよう

---
## 需要モデルの推定と価格弾力性

-   各企業が加入するか否かの意思決定モデル
    $$q(p,x;\theta)=Prob(join|p,x;\theta)=\frac{\exp\left(\alpha(x)+\beta(x)p\right)}{1+\exp\left(\alpha(x)+\beta(x)p\right)}$$

    -   $p$：価格、 $x$: 企業の属性（企業規模、職種などなど）


- $\alpha(x)$ がサービスから得られる効用、 $\beta(x)$ が 価格に対する反応度

- 両パラメタは、企業の属性に依存する
$$\alpha(x) = \sum_{k=1}^K \alpha^k x_k, \ \beta(x) = \sum_{k=1}^K \beta^k x_k,$$
    - 変数 $K$ が大量にある。-> 機械学習
    

    
---
## 推定と価格弾力性

- 推定は、変数選択を考慮した機械学習的手法＋ベイズ推定


-   価格弾力性：
    -   99ドル：-0.33 (95%信頼区間: \[-0.41,-0.26\] )
    -   249ドル：-0.82
    -   399ドル：-1.15


---
## 需要推定に基づく最適価格

-   推定した需要関数を用いて、一時点の収入を最大化する最適価格を計算した。

.center[
<img src = "figuretable/fig_optimal_price.png" width="1000">
]

-   なお、継続率を加味した現在価値最大化の価格は261ドル。

---
## 実験１のその後

-   実験２においては、実験１＋需要推定から得られた価格が「本当に収入増加につながるのか？」について、別途実験を行って検証した。
    -   前ページはin-sampleでの予測。その妥当性について、out-of-sampleで検証


-   その後、Ziprecruiter.comは月249ドルに値上げを行った。
    -   249ドルというのは経営陣の意向を踏まえたもの。



