class: title-slide-section, center, middle
name: logistics
# はじめに

---
## 実証分析例：カリフォルニア州の学校データ

- 問：生徒と教師の比率が試験の成績に与える影響について分析しよう。


- `AER` パッケージに入っているカリフォルニア州の学校データを用いる。
    - 詳細は[このページ](https://www.rdocumentation.org/packages/AER/versions/1.2-6/topics/CASchools)を参照のこと。
    
    
- 以下、２ステップで分析を進めていこう。
    - Step 1：下準備から記述的分析
    - Step 2：回帰分析


- 謝辞：このスライドは [Introduction to Econometrics with R](https://www.econometrics-with-r.org/index.html) に基づく。



---
class: title-slide-section, center, middle
name: logistics
# Step 1：下準備から記述的分析

---
## パッケージの準備

- 以下のパッケージを用いる。
    - `AER`：データソース
    - `dplyr`：データの処理
    - `summarytools`：記述統計の表示
    - `stargazer`：結果表示全般
    - `fixest` ：計量経済学的な推定

    
---

``` {r, warning = FALSE, message = FALSE}
# clear workspace
rm(list =ls())

# load packages
library(AER)
library(tidyverse)
library(summarytools)
library(stargazer)
library(fixest)
```

---
## 注：プログラミングのTips

- Rで外部ライブラリを利用する際には、以下の二通りがある。
  - `library()`で事前に読み込み、パッケージ内の関数を呼び出す。例：require(dplyr)の後にmutate().
  - パッケージ全体を事前に読み込みはせず、特定のパッケージの関数を使う際に[package_name]::[function_name]として呼び出す。例：dplyr::mutate().


- 関数の使う際に、`function(XX, YY)` と `function(input1 = XX, input2 = YY) `という二通りの書き方がある。
  - 前者の場合は、関数で指定されている順に引数を書く必要がある。
  - 後者の場合はその必要はない。例えば、`function(input2 = YY, input1 = XX) ` としてもOK.



---
## データの読み込み

- `AER` パッケージに入っているデータ`CASchools`を読み込む。

```{r}
# ワークスペースにデータを読み込む。
data(CASchools)
```

- データフレーム`CASchools`がEnviornmentの中に登場する。

---
### データの概観

- `head()` 関数でデータの概観を見る。

```{r}
head(CASchools)

```

- `View()` 関数を使って、データセット全体をウィンドウで見ることもできる。

---
## 変数の作成

- `dplyr` パッケージを用いて、必要な変数を作成する。
  - STR：生徒と学生の比率
  - score：試験点数の平均


```{r}
CASchools <- CASchools %>% 
    mutate(STR = students / teachers) %>% 
    mutate(score = (read + math) / 2) 
```

- コメント：パイプ演算子 `%>%` と 代入演算子 `->` に慣れよう。(詳しくは口頭で)

---
## 記述統計

- 記述統計（descriptive statistics）を作る方法はいろいろ。

- `summary()` 関数:Base Rのもの。見た目がイマイチ。

- `stargazer` パッケージ：個人的好み。しかし開発停止によりいろいろ不具合。

- `summarytools` パッケージ：比較的新しいパッケージ

---
## `summary()` 関数


```{r}
summary(CASchools)
```


---
## (参考) 手書きバージョン

```{r}
# STR と score の標本平均の計算
avg_STR <- mean(CASchools$STR) 
avg_score <- mean(CASchools$score)

# STR と score の標本標準偏差の計算
sd_STR <- sd(CASchools$STR) 
sd_score <- sd(CASchools$score)

# 分位点の設定と計算
quantiles <- c(0.10, 0.25, 0.4, 0.5, 0.6, 0.75, 0.9)
quant_STR <- quantile(CASchools$STR, quantiles)
quant_score <- quantile(CASchools$score, quantiles)

# データフレームに集計
DistributionSummary <- data.frame(Average = c(avg_STR, avg_score), 
                                  StandardDeviation = c(sd_STR, sd_score), 
                                  quantile = rbind(quant_STR, quant_score))
```

---
```{r}
head(DistributionSummary)
```


---
## `stargazer` パッケージ
- 個人的なおすすめは `stargazer::stargazer()` 関数を使うこと。

```{r}
stargazer(CASchools, 
          type = "text")
```

---

- 報告したい要約統計を選択できる。

```{r}
stargazer(CASchools,
          type = "text", 
          summary.stat = c("n", "p75", "sd")) 
```

---
## `stargazer` パッケージの補足

- 詳細は[このページ](https://www.jakeruss.com/cheatsheets/stargazer/#the-default-summary-statistics-table) を参照。

- `stargazer` パッケージは回帰分析の結果を報告するのにも使われる。

---
## `summarytools::descr()` 関数

```{r, eval = FALSE}
descr(CASchools, 
      style = "simple",
       transpose = TRUE)
```

---

```{r, echo = FALSE}

descr(CASchools, 
      style = "simple",
       transpose = TRUE)

```

---
## `summarytools::descr()` 関数

- こちらも報告したい要約統計を選択できる。

```{r, warning = FALSE, message = FALSE}
descr(CASchools,
      style = "simple",
      stats = c("n.valid", "q3", "sd"),
      transpose = TRUE)   
```

---
## 散布図

- 生徒と教師の比率と試験の点数の相関を見よう。 

```{r, eval = FALSE}

plot(score ~ STR, 
     data = CASchools,
     main = "Scatterplot of TestScore and STR", 
     xlab = "STR (X)",
     ylab = "Test Score (Y)")

```

---
.center[
```{r, fig.height=8, fig.width=12, echo=FALSE}

plot(score ~ STR, 
     data = CASchools,
     main = "Scatterplot of TestScore and STR", 
     xlab = "STR (X)",
     ylab = "Test Score (Y)")

```
]

---
## 相関係数

- 数値ベクトルの相関を計算するには、 `cor()` 関数を用いる。
```{r}

cor(CASchools$STR, CASchools$score)

``` 

---
class: title-slide-section, center, middle
name: logistics
# Step 2：回帰分析

---
## 回帰分析のパッケージ

- パッケージの例（その他にもいろいろ）
  - `lm()` 関数: Base R 組み込み。標準誤差の計算が均一分散を仮定。
  - `estimatr` パッケージ：不均一分散における標準誤差を考慮 `lm_robust()`など。
  - `lfe` パッケージ： `felm`関数。大量の固定効果の処理が得意。。
  - `fixest` パッケージ：`feols`関数。 `lfe`よりも計算が早い。操作変数法もOK。


- 本講義では最初から`fixest` パッケージを利用する。


---
## 回帰モデル


- 以下の回帰式を考えよう。
$$score_i = \beta_0 + \beta_1 size_i + \epsilon_i$$
ただし、 $size_i$ はクラスサイズ（生徒と教師の比率）である。
    - 以下、変数 `STR` をクラスサイズとよぶ。
    
---
## `feols`による推定

```{r}
# 変数`STR`の名前付け直し
CASchools <- CASchools %>% 
    dplyr::rename(size = STR) 

# モデルの推定と、変数`model1_summary`への結果の格納
model1_summary <- feols(fml = score ~ size, 
                       data = CASchools,
                       vcov = "hetero")
```

- fml: 回帰式の指定。lmと同じ。なお、固定効果や操作変数を考えると複雑になる（来週以降）
- data: 利用するデータフレームの指定。
- vcov: 標準誤差の計算方法の指定。クラスターに頑健な標準誤差でも利用（来週以降）

---
## 推定結果

- `summary`で結果を表示する。
  
```{r}
summary(model1_summary)
```


---
## 結果の解釈

```{r}
summary(model1_summary)
```

- 解釈
    - 教師一人あたりの生徒が1人増えると、試験の点数が2.2点下がる。
    - 標準誤差はとても小さい。クラスサイズが試験の点数に与える影響は統計的に有意である。
    - (自由度を調整していない)$R^2 = 0.051$ であり、被説明変数の分散のうち5.1%がモデルによって説明されると示唆される。
    - 問：影響の経済的な大きさについて考えてみよう。


---
## 標準誤差の比較

- summaryのvcovに"iid"を指定すると、均一分散における標準誤差を報告する。

```{r}
summary(model1_summary, vcov = "iid")
summary(model1_summary, vcov = "hetero")
```

- 標準誤差の値は、均一分散の仮定 ＜ 不均一分散を許す

---
## 説明変数を追加

詳しくは [こちら](https://www.econometrics-with-r.org/7-6-analysis-of-the-test-score-data-set.html) を参照。

```{r}
# 異なるモデルの推定
spec1 <- feols(score ~ size, data = CASchools, vcov = "hetero")
spec2 <- feols(score ~ size + english, data = CASchools, vcov = "hetero")
spec3 <- feols(score ~ size + english + lunch, data = CASchools, vcov = "hetero")
spec4 <- feols(score ~ size + english + calworks, data = CASchools, vcov = "hetero")
spec5 <- feols(score ~ size + english + lunch + calworks, data = CASchools, vcov = "hetero")
```

---
## 推定結果の表の作成


- **【重要】** summary() で出した結果は自分で見る分にはOKだが、レポートするものとしては不十分。


- 宿題やレポート、論文において、決してsummary()の結果を貼り付けないこと！


- 推定結果は表としてレポートするのが良い。


- 推定結果を表でまとめる方法もいろいろある。
  - `stargazer`: 最近のパッケージ(estimatr::lm_robustやfixest:feols)を受け付けない。
  - `modelsummary`: フレキシブル。その分慣れが必要
  - `fixest::etable`: fixestパッケージに付属するもの。
  
  
- 以下では`fixest::etable`を用いる。

---
## 推定結果をまとめた表を作成

```{r}
result <- etable( list(spec1, spec2, spec3, spec4, spec5),  
                se = "hetero",
                fitstat = c("r2", "n" ) , 
                signif.code = NA, 
                digits = 2, digits.stats = 2,
                depvar = FALSE)

```

- 注：`dict = c(size = "教員・学生比率", `(Intercept)` = "定数項")`オプションで変数名を日本語に変更可能。
- `signif.code =NA` で統計的有意性の記号を外している。昨今の論文は記号を載せない傾向。

---
## 推定結果

```{r, eval = TRUE}
print(result)
```

---
## 練習問題


- 問い：説明変数を増やすとクラスサイズの係数が減少するのはなぜだろうか？（ヒント：欠落変数バイアス）
