Rでのパネルデータ加工

# はじめに
久保研介研究会のサブゼミの資料です！一応、限定公開にしてます。（2020.08. 一般公開しました。）
Rでのパネルデータ加工が資料として纏まっていないので、オリジナルで作成しました。
雑なコードが多いですが、日々Updateしていきます。わかりにくい説明等あれば、どしどし連絡ください！
※基本的な進め方はお兄さん（ゼミの先輩）の資料を参考に作成しました！ありがとうございます！感謝です！

### お題の確認
今回のゴールは経済主体が特許取得者，観察期間が2010年〜2016年のパネルデータ（df_final）を作成することです！元ネタは**[三田際論文製薬班'19](https://43c49aa2-d9d5-427b-b96b-9dc3dd5a6113.filesusr.com/ugd/a2f719_131babca8e0e4d878e31b79b62fbec80.pdf)**で使用したデータです。（クリックすると久保研介研究会のHPへ飛びます。）データは[こちら](https://drive.google.com/drive/folders/1cCFYE0q98kq-UDvwsTqjQhNCeTR2v9GG?usp=sharing)からです。悪用はしないでください。
### 変数の説明
- 特許取得者
 - 特許を取得した人（このデータの経済主体）
- 年度
 - 観察期間（2010年~2016年）
- 特許取得数
 - 当該年に特許を取得した数
- 会社名
 - 特許取得社は在籍している会社名
- 有形固定資産〜研究開発集約度合
 - 在籍している会社の財務情報
- 合併
 - 当該年度において在籍している会社が合併を実施した回数

<img width="666" alt="スクリーンショット 2020-06-27 21.58.13.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/360143/0e4467c5-9555-0bfa-f0a7-87e94ee244b7.png">

### パネルデータ作成のTips
**ゴールまでの道のりを解像度を高く描くことがポイント．**
全体像を描いてから，より細かいデータの動かし方を考えるという流れです．
最初のうちは，いきなりRと向き合うのではなく，**絵をかいたり，言葉に起こしたり**すると良いです．


### 今回の大まかな流れ
Step1: とりあえずデータをくっつけていく!

- 1-1: df_patent_allにdf_patent_medicを結合する：誰がいつどの企業で特許を取得したのか？
- 1-2: df_patent_2にdf_finを結合する：企業の財務状況を追加
- 1-3: df_patent_3にdf_maを結合する：m&a情報を追加する

Step2: 観察主体が特許取得者、観察期間が2010年〜2016年のパネルデータにしていく

- 2-1: パネルデータ化
- 2-2: パネルデータにしたことで生まれた欠損値を処理していく

# 使うデータの読み込み
では早速始めていきましょう！
以下の4つのデータを使ってパネルデータを作成していきます．
xlsxファイルをインポートする際にopenxlsxというライブラリを使用しています．ダウンロードが出来ていない人がいれば，```install.packages('openxlsx')```と入力してください．

```r
library(openxlsx)
df_patent_medic <- read.xlsx("製薬論文データ　製薬企業の特許一覧.xlsx")
df_patent_all <- read.table("inventor.txt", header = T,sep = ",")
df_ma <- read.xlsx("製薬論文データ　M&A一覧.xlsx")
df_fin <- read.xlsx("製薬論文データ　企業財務データ.xlsx")
```

### 変数の説明（わかりにくいものだけに留めます）
**df_patent_all**

- ida 
 - 特許開発番号
- seq
 - 複数の開発者がいたことを示す
- ida_seq
 - idaにseqの値を加えたもの

**df_fin**

- asset
 - 有形固定資産額
- ope
 - 営業利益
- sale
 - 売上高
- con
 - 研究開発集約度

# Step1: とりあえずデータをくっつけていく!
## Step1-1: df_patent_allにdf_patent_medicを結合する：誰がいつどの企業で特許を取得したのか？
両方のデータに共通する特許開発番号で縦方向に結合してみます。
### データの型の確認
データを加工する際には、データの型を把握しておく必要があります。
データの型によって、処理方法が違うためです。例えば、見た目が数字であっても文字列と認識されていれば、平均値などの計算をすることができません。

```r
class(df_patent_medic$企業名)
sapply(df_patent_medic, class)
sapply(df_patent_all, class)
```
#### [データの型の種類](https://stats.biopapyrus.jp/r/basic/data-type.html)※クリックすると参考サイトに飛びます。
- numeric
 - 数値
- integer
 - 整数
- character
 - 文字列
- factor
 - 因子（文字列に順番が与えられたもの）

### 問題点：特許番号がちょいと違うのでdf_patentのidaという変数のデータに揃える
結合するためには、同じ値である必要があるので、ハイフンが邪魔ですね。ここでは**文字列の置換**を行うための **[sub](https://stats.biopapyrus.jp/r/basic/string.html)**というコマンドを使用して、処理しています。
また、文字列ではなく数値として扱いたいので、 **[as.character](https://stats.biopapyrus.jp/r/basic/data-type.html)**で数値に変換しています。**as.変換したいデータの型**というイメージです。

```r
df_patent_medic$特許出願番号<- sub("-", "", df_patent_medic$特許出願番号)
df_patent_all$ida <- as.character(df_patent_all$ida)
class(df_patent_all$ida)
```
### データの結合
データの結合をするためには、 **[merge](http://cse.naro.affrc.go.jp/takezawa/r-tips/r/43.html)**というコマンドを使用します。2行目のコードでは、必要な列だけを **[抽出](https://qiita.com/gigatune/items/f3aa0afef7f50ab791cd)**しています。**[行の指定,列の指定]**というイメージです。

```r
df_patent <- merge(df_patent_all, df_patent_medic, by.x = "ida", by.y = "特許出願番号")
df_patent <- df_patent[, c("name", "年度", "企業名")]
```
#### データ結合のTips
データを結合する際に**気をつけておきたいポイント**を紹介しておきます。

- データをどの方向で結合するのか？
 - 縦方向ならば、**[rbind](https://webbeginner.hatenablog.com/entry/2015/02/06/132256)**を、横方向ならば、**[merge](http://cse.naro.affrc.go.jp/takezawa/r-tips/r/43.html)**を使用します。
- 横方向で結合する場合には、それぞれのデータの**どの変数をkey**として結合するのか？
 - keyとする変数が同じ値であるか？（データの**見た目**、**データの型**）
 
### データを見て考えること①：1年に複数の特許を取得している場合
こういう場合はデータ上ではどのように現れるんだろう？というような思考ができれば、データ加工スキルが上がると思います。
例えば、1年に複数の特許を取得しているケースは、どのようにデータ上で表現されるでしょうか？おそらく、複数行で表現されているはずです。これは、のちのちパネルデータにしていく際に困りますね。（パネルデータでは、各観察期間につき各経済主体のデータは1つずつです。）
そのため、まずは各特許開発者が各年度において、いくつ特許を取得しているかを集計したいと思います。データの集計は **[dplyr](https://qiita.com/matsuou1/items/e995da273e3108e2338e)**の```group_by```と```summarise```を組み合わせることで出来ます。dplyrはデータ加工のためのライブラリで非常に便利です。各自URLから詳しい使い方をチェックしてみてください。

```r
install.packages("dplyr")
library(dplyr)
df_patent_1 <- df_patent %>% 
  group_by(name, 年度, 企業名) %>% 
  summarise(patent = n())
```
はじめは慣れない書き方だと思いますが、df_patentとdf_patent_1のデータを見比べて、どのような変化がおこったのか？を確認してみてください！ちなみに ```%>%``` (パイプライン)については、この **[サイト](https://www.medi-08-data-06.work/entry/paipe_r)**を参考にしてください。

### データを見て考えること②：会社を移籍している場合
もう一つ、特許開発者が観察期間内で移籍している場合も考えられます。（寒川　賢治さんとか）移籍している人は、同じ年度に企業名が複数あることが考えられます。そのため、特許開発者と年度を軸に行数を集計してみると移籍しているかどうかがわかると思います。

```r
df_check <- df_patent_1 %>% 
  group_by(name, 年度) %>% 
  summarise(移籍 = n())
```
今回は会社を移籍しているケースは考慮しないこととします。なので、特許開発者と年度に重複がある場合は重複しているどちらかを削除します。**データの重複削除**は **[distinct](https://a-habakiri.hateblo.jp/entry/2016/11/29/215013)**というコマンドを使うと出来ます。

```r
df_patent_2 <- df_patent_1 %>% 
  arrange(name, 年度, 企業名) %>% 
  distinct(name, 年度, .keep_all = T)
```

## Step1-2: df_patent_2にdf_finを結合する：企業の財務状況を追加
### データの型の確認
Step1-1と同様に、df_finというデータを結合していきましょう。今回は**企業名と年度が結合のkey**になります。結合に入る前に、まずはデータの型を確認しておきましょう。

```r
sapply(df_fin, class)
```

### 問題点：会社名がちょいと違うのでdf_fin側に揃える
結合のkeyである企業名が両方のデータで共通しているかどうかを確認してみましょう。```unique```というコマンドで名前通り**ユニークな値**を取得することが出来ます。

```r
unique(df_patent_2$企業名)
unique(df_fin$company)
```
田辺三菱製薬株式会社、アステラス、大日本住友の3社の名前が異なりますね。今回はdf_finでの企業名にそろえましょう。正しく揃えられたかどうかを最後のコードで確認しています。```%in%```の左の値が右に含まれていたらTRUEを返してくれるコードになります。実行すると、全てTRUEと表示されるので、正しく企業名を揃えられたことが確認出来ます。

```r
df_patent_2$企業名 <- sub("田辺三菱製薬株式会社", "田辺三菱製薬", df_patent_2$企業名)
df_patent_2$企業名 <- sub( "アステラス", "アステラス製薬" , df_patent_2$企業名)
df_patent_2$企業名 <- sub("大日本住友", "大日本住友製薬", df_patent_2$企業名)
unique(df_patent_2$企業名) %in% unique(df_fin$company)
```

### データの結合
keyとなる変数を揃えることが出来たので、データの結合に取り掛かりましょう。今回は先ほどと違い、**企業名と年度の2つの変数がkey**となります。そのため、```by.x```と```by.y```での指定をベクトル表記の```c()```を用いて複数指定します。結合したデータをdf_patent_3とします。

```r
df_patent_3<- merge(df_patent_2, df_fin, by.x = c("年度",  "企業名"), by.y = c("X2", "company") )
```

## Step1-3: df_patent_3にdf_maを結合する：m&a情報を追加する
### 同様に企業名のチェック
最後のデータ結合になります。前回と同様に**企業名と年度が結合のkey**となります。だんだん要領を掴んできたのではないでしょうか？いつも通り各自でデータの型を確認しておいてください。企業名が揃っているかどうかのチェックもしておきましょう。

```r
unique(df_ma$企業名) %in% unique(df_patent_3$企業名)
```
### データの集計：先にM＆A件数をカウントしておく
手元のデータでは各企業の同じ年のM&Aは複数行に渡っています。今回、結合したいデータは**当該年に企業が何件M＆Aをしたのかどうか**なので集計してみましょう。復習ですが、**データの集計**は **dplyr**の```group_by```と```summarise```を組み合わせることで出来ましたね。

```r
df_ma_1 <- df_ma %>% 
  group_by(企業名,year) %>% 
  summarise(MAcount = n())
```

### データの結合
さて、今回最後の結合です。先ほどと基本的には同じですが、今回は**x側のデータ（つまりdf_patent_3）を結合後も残しておきたいので、```all.x =TRUE```という指定が必要**になります。この指定をしないで結合をしてしまうと、M&Aの記録がない年のデータが消えてしまいます。実際に試してみると、納得できるはずです。

```r
df_patent_4 <- merge(df_patent_3, df_ma_1, by.x = c("企業名", "年度"),by.y = c("企業名","year"), all.x =TRUE)
```

### 欠損値処理
MAcountがNAとなっている部分はM&Aをしていないということなので、0を代入しておきます。欠損値への代入はこのように書きましたね。

```r
df_patent_4$MAcount[is.na(df_patent_4$MAcount)] <- 0
```

### データの整理
列番号を並び替えることで、列の並び替えもできます。経済主体と観察期間を左にもっていき、パネルデータを意識した形にしておきましょう。

```r
df_patent_4 <- df_patent_4[, c(3, 2, 1, 4, 5, 6, 7, 8, 9)]
```

# Step2: 観察主体が特許取得者、観察期間が2010年〜2016年のパネルデータにしていく

## Step2-1: パネルデータ化
いよいよ終盤ですね。 **[complete](https://notchained.hatenablog.com/entry/2017/04/04/220236)**というコマンドを用いるとパネルデータにすることが出来ます。コードの書き方としては、dataを指定した後に、経済主体を指定。その後に、観察期間を指定します。下のコードのように期間を手動で指定することも可能です。

```r
library(tidyverse)
df_panel <- complete(data = df_patent_4, name, 年度 = 2010:2016)
```

## Step2-2: パネルデータにしたことで生まれた欠損値を処理していく
### どうやって埋めるか？
先に**企業名**を埋めます。企業名が埋まれば、**財務情報もM＆A情報**も埋まるはずです。
でもその前にpatentを埋めちゃいましょう。patentがNAということは特許を開発していないということので、0を代入してあげます。

```r
df_panel$patent[is.na(df_panel$patent)] <- 0
```

### 企業名を数値化
如何せん文字列だと扱いづらいので，企業名を数字に変換します。文字列だと扱いにくい理由は後々わかるはずです。
やりたいことはこんなイメージです。

```r
df_panel$企業名 <- sub("田辺三菱製薬",1, df_panel$企業名)
df_panel$企業名 <- sub("アステラス製薬" ,2, df_panel$企業名)
df_panel$企業名 <- sub("大日本住友製薬", 3, df_panel$企業名)
df_panel$企業名 <- sub("エーザイ",4, df_panel$企業名)
df_panel$企業名 <- sub("協和発酵キリン" ,5,  df_panel$企業名)
df_panel$企業名 <- sub("第一三共", 6, df_panel$企業名)
df_panel$企業名 <- sub("武田薬品工業", 7, df_panel$企業名)

df_panel$企業名 <- as.numeric(df_panel$企業名)
```
これはだいぶめんどくさいですね。これを1行でやると、下ののコードになるかと思います。一度**factor(因子)**に変換することで、**文字列に順番がつきます**。この性質を利用して、文字列を数値(ここでは整数)に変換します。

```r
df_panel$企業名 <- as.integer(as.factor(df_panel$企業名))
```

### 企業名を埋める
では企業名（数字）を埋めていきます。今回は簡略化のため移籍を考慮しないので、観察期間において経済主体は同一の企業に勤めているとします。そのため、**各特許開発者ごとの企業名の最大値（最小値でも良い）を代入してあげれば、欠損値が上手く埋まります。**この手法はパネルデータ作成の特有のやり方だと思います。なかなかすぐに腹落ちしないと思います。私もそうでした。

```r
df_panel<- df_panel %>% 
  group_by(name) %>% 
  mutate(会社名 = max(企業名, na.rm = T))
```
今回は**dplyr**パッケージに入っている**group_by**と**mutate**というコマンドを組み合わせて、集計した値（特許開発者ごとの最大値）を新しい変数（会社名）として結合させました。**mutateを使うと新しい変数を作成する**ことが出来ます。また、```max```というコマンドで最大値を算出しています。**```na.rm = T```と指定することで、欠損値によるエラーを回避できます。**新しく**会社名**という列が出来ていると思うので、確認してみてください。企業番号で埋まっているはずです。
### 財務情報・MAも埋めていく
同様に財務情報とM&A情報も埋めていきます。これらの情報は**企業名と年度ごとに**埋めたいので、先ほど埋めた**会社名と年度**で集計しています。下のコードを実行してみて、欠損値がなくなっているか確認してみましょう。

```r
df_panel<- df_panel %>% 
  group_by(会社名, 年度) %>% 
  mutate(有形固定資産 = max(asset, na.rm = T) , 営業利益 = max(ope, na.rm = T),売上 =max(sale, na.rm = T)
               , 研究開発集約度 = max(con, na.rm = T), 合併 = max(MAcount, na.rm = T))
```

### 最終的なデータに揃える 
必要な変数だけを残します。同じく**dplyr**の```select```というコマンドで、**列を抽出**することが出来ます。また、せっかくなので英語の変数名を日本語にしておきましょう。同じく**dplyr**の```rename```というコマンドで、**列名を変更**することが出来ます。

```r
df_final <- df_panel %>% 
  select(name, 年度, patent,会社名, 有形固定資産,営業利益,売上, 研究開発集約度,合併 ) %>% 
  rename(特許取得者 = name, 特許取得数 = patent)
```
**これで完成です！**このデータを使った分析の詳細については、**[三田際論文製薬班'19](https://43c49aa2-d9d5-427b-b96b-9dc3dd5a6113.filesusr.com/ugd/a2f719_131babca8e0e4d878e31b79b62fbec80.pdf)**を見てください。（クリックすると久保研介研究会のHPへ飛びます。）**DID分析**を使って、**合併が研究開発に与える影響**を分析している論文になります。

# 最後に
データ加工は結構疲れますが、データ分析で一番重要な段階です！Garbage in Garbage out の精神で、コツコツ頑張っていきましょう👍