GES 001 演習8

2024年2月13日

準備

（自分のPCまたは教室のPCに）ログイン
ウェッブ・ブラウザー（Google Chrome など）を起動
- Moodle の GES001 経済と経済学のサイトから、このスライドのページを表示（リンク「Rでデータサイエンス」https://ds-sl.github.io/intro2r/ges001/）
（別のタブまたはウィンドウで）PositCloud にログイン[Posit.cloud]
- アカウントのない人はサイン・アップ [共有プロジェクト] から、Save a Permanent Copy）
- RStudio を自分のコンピュータにインストールしている人は起動
リンクの右上の Raw ボタンの右の Copy a raw file からコピーして演習用 R Markdown ファイルを作成（あとで再度解説します）[Rmd]

ファイル

ファイルを作成するときの注意

Posit Cloud のときは、まず、Login し、intro2rj のプロジェクトに入り、ファイルから、ges001 のフォルダーを選択して、移動します。この中に、data フォルダが作成されていることを確認し、このges001 に、新しく作成したファイルを保存します。新しく作成したファイルの入っているフォルダーの中に、data フォルダがあることが大切です。そこに、データを書き込みます。
RStudio の場合には、まずは、あたらしい　Project を作成します。File > New Project から作成します。すでに、作成してある場合は、それを、Open Project や、Recent Project から開きます。その中に、新しいファイルを作成します。作成したフォルダーに、data フォルダがあることを確認してください。新しく作成したファイルの入っているフォルダーの中に、data フォルダがあることが大切です。そこに、データを書き込みます。

ファイルを提出するときの注意

RStudio の場合には、自分の PC に作成したファイルがありますから、問題ないと思いますが、Posit Cloud で作成した場合には、提出したいファイルの左にあるチェックボックスをチェックします。Files の右端にある、ギアマークの Export を押すと、ダウンロードできます。それを提出してください。末尾が、nb.html となっているものを提出していただくのがよいですが、よくわからないときは、nb.html ファイルと、Rmd ファイルと両方提出してください、

課題1について

課題1：解説リンク
課題1：解答例リンク
探索的データ分析（課題2の解説にかえて）in progress リンク

課題1・課題2 追加提出ボックスについて

期限：2023年2月17日23時59分（特別な理由がある場合は連絡してください）

第8週

02/08(TH)　紛争と貧困の連鎖1

紛争と貧困の連鎖2

第8週と第9週では、戦争、気候変動、COVID-19の複合的要因が世界的な食糧危機を引き起こし、格差と貧困を加速させている現状について、議論します。

02/13(TU)　Rでデータサイエンス8：難民、軍事支出　 [Main]・[授業]

Global Link-Refugees

Refugee population by country or territory of asylum：SM.POP.REFG [Link]

Refugee population by country or territory of origin：SM.POP.REFG.OR [Link]

Net ODA received (% of GNI)：DT.ODA.ODAT.GN.ZS [Link]

Net official development assistance and official aid received (current US$)：DT.ODA.ALLD.CD [Link]

Net ODA received (% of central government expense)：DT.ODA.ODAT.XP.ZS [Link]

Military expenditure (current USD)：MS.MIL.XPND.CD [Link]

Military expenditure (% of general government expenditure)：MS.MIL.XPND.ZS [Link]

Arms imports (SIPRI trend indicator values)：MS.MIL.MPRT.KD [Link]

Arms exports (SIPRI trend indicator values)：MS.MIL.XPRT.KD [Link]

探索的データ分析（EDA）ワークフロー

表題 Title
概要 Abstract
データ情報：データ名、データコード、変数名、データ概要
各項目に、気づいたこと、疑問などを記録（EDAの核となる部分です）
データの取得 - ダウンロードして保存しておく
- 準備：library(tidyverse); library(WDI)
- df_w6eda <- WDI(indicator = c(co2pcap = "EN.ATM.CO2E.PC", forest = "AG.LND.FRST.ZS"), extra = TRUE)
- 二つ以上の指標、それぞれに変数名 co2pcap, forest
- 追加情報：extra = TRUE; region 地域情報, income 収入レベル
- データの確認：データが十分あるか、国名リスト、地域名リスト
分析する国のリスト：南部アフリカ関税同盟、国や地域を選択

探索的データ分析（Exploratory Data Analysis）ワークフロー（続き）

分析・視覚化
- 各年毎のデータの数の棒グラフ
- 経年変化
  - 日本（身近な馴染みのあるデータを確認）
  - 南部アフリカ関税同盟、選択した国・地域
- 分布
  - ヒストグラム（年を指定した度数分布）
  - 値が大きい10カ国とその値の棒グラフ
  - 値が小さい10カ国とその値の棒グラフ
- 相関
  - 散布図（scatter plot）x-軸と、y-軸に、それぞれの変数をとる
  - 相関係数

演習 2月13日（火）

変形（transformation）

行の選択：filter(country %in% c("Japan", "China"), distinct(country), drop_na(value)
列の選択：select(country, iso2c, year, var1, var2, region)
横並びの変数を縦並びに：name 列と、value 列を作り、var_1, var_2, … , var_n を name 列に、対応する値を　value 列に
- pivot_longer(cols = c(var_1, var_2, ... , var_n))
- pivot_longer(cols = c(female_unemploy, male_unemploy))
- pivot_longer(cols = c(electricity_urban, electricity_rural))
- pivot_longer(cols = electricity_urban:sanitation_rural))
活用：pivot_longer(cols = c(var1, var2))|> ggplot(aes(year, value, col = name) + geom_line()

箱ひげ図（boxplot）

カテゴリ変数（country, region, income など）毎の分布を描く

ggplot(aes(categorical_var, numerical_var)) + geom_boxplot()
ggplot(aes(region, electricity_urban)) + geom_boxplot()
ggplot(aes(income, electricity_rural)) + geom_boxplot()
pivot_longer(cols = c(var1, var2)) |> ggplot(aes(name, value)) + geom_boxplot()
year は、整数なので数値変数（numerical variable）。これをカテゴリー変数として扱うためには、factor(year) を使う。x 軸のラベルなどはあとから修正が必要。
- filter(year %in% c(1960,1990,2020))|> ggplot(aes(factor(year), var)) + geom_boxplot()

内容

Refugee population by country or territory of asylum：SM.POP.REFG [Link]
Refugee population by country or territory of origin：SM.POP.REFG.OR [Link]
Net ODA received (% of GNI)：DT.ODA.ODAT.GN.ZS [Link]
Net official development assistance and official aid received (current US$)：DT.ODA.ALLD.CD [Link]
Net ODA received (% of central government expense)：DT.ODA.ODAT.XP.ZS [Link]
Military expenditure (current USD)：MS.MIL.XPND.CD [Link]
Military expenditure (% of general government expenditure)：MS.MIL.XPND.ZS [Link]
Arms imports (SIPRI trend indicator values)：MS.MIL.MPRT.KD [Link]
Arms exports (SIPRI trend indicator values)：MS.MIL.XPRT.KD [Link]

データの取得

準備

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(WDI)

df_peace <- WDI(
  indicator = c(refugee_asylum = "SM.POP.REFG",
                refugee_origin = "SM.POP.REFG.OR",
                oda_gni = "DT.ODA.ODAT.GN.ZS",
                oda_usd = "DT.ODA.ALLD.CD",
                oda_gov = "DT.ODA.ODAT.XP.ZS",
                military_usd = "MS.MIL.XPND.CD",
                military_gov = "MS.MIL.XPND.ZS",
                arms_imports = "MS.MIL.MPRT.KD",
                arms_exports = "MS.MIL.XPRT.KD"), extra = TRUE)

write_csv(df_peace, "data/peace.csv")

df_peace <- read_csv("data/peace.csv")

## Rows: 16758 Columns: 21
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr   (7): country, iso2c, iso3c, region, capital, income, lending
## dbl  (12): year, refugee_asylum, refugee_origin, oda_gni, oda_usd, oda_gov, ...
## lgl   (1): status
## date  (1): lastupdated
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

データの確認

df_peace

データの変数（列）の確認

str(df_peace)

## spc_tbl_ [16,758 × 21] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ country       : chr [1:16758] "Afghanistan" "Afghanistan" "Afghanistan" "Afghanistan" ...
##  $ iso2c         : chr [1:16758] "AF" "AF" "AF" "AF" ...
##  $ iso3c         : chr [1:16758] "AFG" "AFG" "AFG" "AFG" ...
##  $ year          : num [1:16758] 2014 1971 2006 2013 1995 ...
##  $ status        : logi [1:16758] NA NA NA NA NA NA ...
##  $ lastupdated   : Date[1:16758], format: "2023-12-18" "2023-12-18" ...
##  $ refugee_asylum: num [1:16758] 300421 NA 34 16861 19605 ...
##  $ refugee_origin: num [1:16758] 2596259 NA 2107510 2556483 2679132 ...
##  $ oda_gni       : num [1:16758] 24.04 2.38 NA 25.5 NA ...
##  $ oda_usd       : num [1:16758] 4.94e+09 4.44e+07 2.90e+09 5.15e+09 2.13e+08 ...
##  $ oda_gov       : num [1:16758] 54.1 NA 201.8 60.8 NA ...
##  $ military_usd  : num [1:16758] 2.68e+08 NA 1.31e+08 2.17e+08 NA ...
##  $ military_gov  : num [1:16758] 5.11 NA 10.36 4.31 NA ...
##  $ arms_imports  : num [1:16758] 3.17e+08 1.49e+08 3.00e+06 2.55e+08 NA 3.30e+07 NA NA NA NA ...
##  $ arms_exports  : num [1:16758] NA NA NA NA NA NA NA NA NA NA ...
##  $ region        : chr [1:16758] "South Asia" "South Asia" "South Asia" "South Asia" ...
##  $ capital       : chr [1:16758] "Kabul" "Kabul" "Kabul" "Kabul" ...
##  $ longitude     : num [1:16758] 69.2 69.2 69.2 69.2 69.2 ...
##  $ latitude      : num [1:16758] 34.5 34.5 34.5 34.5 34.5 ...
##  $ income        : chr [1:16758] "Low income" "Low income" "Low income" "Low income" ...
##  $ lending       : chr [1:16758] "IDA" "IDA" "IDA" "IDA" ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   country = col_character(),
##   ..   iso2c = col_character(),
##   ..   iso3c = col_character(),
##   ..   year = col_double(),
##   ..   status = col_logical(),
##   ..   lastupdated = col_date(format = ""),
##   ..   refugee_asylum = col_double(),
##   ..   refugee_origin = col_double(),
##   ..   oda_gni = col_double(),
##   ..   oda_usd = col_double(),
##   ..   oda_gov = col_double(),
##   ..   military_usd = col_double(),
##   ..   military_gov = col_double(),
##   ..   arms_imports = col_double(),
##   ..   arms_exports = col_double(),
##   ..   region = col_character(),
##   ..   capital = col_character(),
##   ..   longitude = col_double(),
##   ..   latitude = col_double(),
##   ..   income = col_character(),
##   ..   lending = col_character()
##   .. )
##  - attr(*, "problems")=<externalptr>

変形

変数の選択（selecting）

df_pc <- df_peace |> 
  select(country, iso2c, year, refugee_asylum:arms_exports, region, income)
df_pc

縦長の表（Long Table）に変換

df_pc_long <- df_pc |> 
  pivot_longer(refugee_asylum:arms_exports)
df_pc_long

各年毎のデータの数の確認

df_pc_long |> drop_na(value) |> 
  ggplot(aes(year)) + geom_bar()

各指標ごとに積み上げ棒グラフにしてみます

fill を col に変えるとどうなるでしょうか。

df_pc_long |> drop_na(value) |> 
  ggplot(aes(year, fill = name)) + geom_bar()

気づいたこと、疑問

国と地域

地域の `iso2c` コード

country には、国と地域両方が入っています。地域の iso2c は以下のものです。

REGION <- c("1A", "1W", "4E", "6F", "6N", "6X", "7E", "8S", "A4", "A5", 
"A9", "B1", "B2", "B3", "B4", "B6", "B7", "B8", "C4", "C5", "C6", 
"C7", "C8", "C9", "D2", "D3", "D4", "D5", "D6", "D7", "EU", "F1", 
"F6", "M1", "M2", "N6", "OE", "R6", "S1", "S2", "S3", "S4", "T2", 
"T3", "T4", "T5", "T6", "T7", "V1", "V2", "V3", "V4", "XC", "XD", 
"XE", "XF", "XG", "XH", "XI", "XJ", "XL", "XM", "XN", "XO", "XP", 
"XQ", "XT", "XU", "XY", "Z4", "Z7", "ZB", "ZF", "ZG", "ZH", "ZI", 
"ZJ", "ZQ", "ZT")

地域のリストを表示

df_pc |> filter(iso2c %in% REGION) |> distinct(country, iso2c)

国名とその地域・所得レベルを表示

df_pc |> filter(!(iso2c %in% REGION)) |> distinct(country, iso2c,region,income)

国や地域の選択

ASEAN <- c("BN", "KH", "ID", "LA", "MY", "MM", "PH", "SG", "VN")

BRICS <- c("Brazil", "Russian Federation", "India", "China", "South Africa")

SACU <- c("South Africa", "Namibia", "Eswatini", "Botswana", "Lesotho")

CHOSEN_COUNTRIES <- c("Suriname", "Belize", "Brazil", "Colombia")

変数同士の相関関係: NA ではないところのみ選択

df_peace |> filter(!(iso2c %in% REGION)) |> # 地域以外（国のみ選択）
  drop_na(refugee_asylum:arms_exports) |> 
  select(refugee_asylum:arms_exports) |> cor()

##                refugee_asylum refugee_origin      oda_gni     oda_usd
## refugee_asylum     1.00000000     0.09529925  0.277803615  0.19080685
## refugee_origin     0.09529925     1.00000000 -0.055259672  0.02793923
## oda_gni            0.27780361    -0.05525967  1.000000000  0.18238369
## oda_usd            0.19080685     0.02793923  0.182383686  1.00000000
## oda_gov            0.20651969    -0.06003550  0.934039951  0.15519192
## military_usd      -0.04219972    -0.03483676 -0.297119009  0.45219438
## military_gov       0.10960070    -0.08599254 -0.002842719 -0.19992294
## arms_imports      -0.08180642    -0.11327444 -0.208204700  0.46577548
## arms_exports       0.15657351     0.28715525 -0.066830867  0.16262742
##                    oda_gov military_usd military_gov arms_imports arms_exports
## refugee_asylum  0.20651969  -0.04219972  0.109600697 -0.081806419   0.15657351
## refugee_origin -0.06003550  -0.03483676 -0.085992543 -0.113274438   0.28715525
## oda_gni         0.93403995  -0.29711901 -0.002842719 -0.208204700  -0.06683087
## oda_usd         0.15519192   0.45219438 -0.199922938  0.465775477   0.16262742
## oda_gov         1.00000000  -0.27186566 -0.015283892 -0.174893747  -0.09637824
## military_usd   -0.27186566   1.00000000 -0.186778574  0.659515539   0.02278056
## military_gov   -0.01528389  -0.18677857  1.000000000 -0.005298621  -0.01080562
## arms_imports   -0.17489375   0.65951554 -0.005298621  1.000000000  -0.13894503
## arms_exports   -0.09637824   0.02278056 -0.010805616 -0.138945034   1.00000000

値の正負：正の相関・負の相関、1に近いと強い正の相関、-1に近いと強い負の相関

相関係数：散布図の近似（回帰）直線の傾きが正なら正、負なら負、直線に近い程、1 または-1 に近い

経年変化

a. 日本

df_pc_long |> filter(country == "Japan") |> drop_na(value) |>
  ggplot(aes(year, value, col = name)) + geom_line()

df_pc_long |> filter(country == "Japan") |> drop_na(value)

df_pc_long |> filter(iso2c %in% ASEAN) |> drop_na(value) |>
  ggplot(aes(year, value, col = country, linetype = name)) + geom_line()

df_pc_long |> filter(iso2c %in% ASEAN) |> drop_na(value) |> 
  filter(value < 80) |>
  ggplot(aes(year, value, col = country, linetype = name)) + geom_line()

df_pc_long |> filter(country %in% SACU) |> drop_na(value) |>
  ggplot(aes(year, value, col = country, linetype = name)) + geom_line()

箱ひげ図（Boxplot）

df_pc_long |> filter(name == "refugee_asylum") |> 
  filter(region != "Aggregates") |> drop_na(value) |>
  ggplot(aes(region, value)) + geom_boxplot() + coord_flip()

df_pc_long |> filter(name == "refugee_asylum") |> 
  filter(region != "Aggregates") |> drop_na(value) |>
  ggplot(aes(region, value)) + geom_boxplot() + scale_y_log10() + coord_flip()

df_peace |> filter(region != "Aggregates") |> drop_na(refugee_asylum) |>
  ggplot(aes(region, refugee_asylum)) + geom_boxplot() + scale_y_log10() + coord_flip()

df_pc_long |> filter(name %in% c("refugee_asylum", "refugee_origin")) |> 
  filter(income != "Aggregates") |> drop_na(value) |>
  ggplot(aes(income, value, fill = name)) + geom_boxplot() + scale_y_log10() + coord_flip()

df_pc_long |> filter(income != "Aggregates", income != "Not classified") |> drop_na(value) |>
  ggplot(aes(income, value, fill = name)) + geom_boxplot() + scale_y_log10() + coord_flip()

## Warning in self$trans$transform(x): NaNs produced

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 835 rows containing non-finite values (`stat_boxplot()`).

ファイルリンク

基本的には、PositCloud（https://posit.cloud/）を使って実習

探索的データ分析（EDA） -
- 課題1+2:例平和（peace.Rmd） [リンク], [Rmd]
- 課題1+2：テンプレート：w8eda_template.Rmd [リンク], [Rmd]

2023.2.17.23:59 までに Moodle の課題２提出ボックスに提出してください。

できたグラフだけでも構いませんから、期限までに提出してください。解説を加え、課題２再提出ボックスを作成する予定です。

課題2

以下の指標の中から、二つ以上（複数）を選択して、データの概要（description）を記録し、データを WDI で取得し、以下の分析をする。

1 各年毎のデータの数の棒グラフ

2 経年変化を表す折れ線グラフ a. 日本、またはデータがある国 b. 選択したいくつかの国

3 複数の指標の値を一列に含む縦長の表（Long Table）を作成し　 a. 日本、またはデータがある国の、複数の指標を色分けした経年変化のグラフ　 b. 選択したいくつかの国についての経年変化のグラフを、国を色分けし、指標は線の種類を変えたグラフ

課題2 つづき

4 二つのデータの散布図- NA は取り除くこと。（log10 スケールを用いる場合は値が正のもののみに限定） a. （地域を除き）国のみの散布図（近似（回帰）直線を表示） b. 最近の年の（地域を除き）国のみの散布図　（近似（回帰）直線を表示） c. b に対応する相関係数

5 カテゴリー変数（Categorical Variable: region, income, year など）と、数値変数（Numberical Variable）一組についての箱ヒゲ図（Boxplot）

それぞれについて考察（気づいたこと、疑問など）を記す

2023.2.17.23:59 までに Moodle の課題２追加提出ボックスに提出してください。

確認作業

Preview で確認。
Web Browser で、w5_c123456.nb.html など、R Notebook を見て確認。
もし、問題があれば、Run ボタンの右の三角から、Run All を選択し、エラーがでないか確認。
最初にもどる。

エラー対応

入力したときには、例を参照して、スペルを確認。全角や変数名の空白は問題を起こす。括弧　() や引用符がペアでマッチしているか、確認。引用符が入っていなかったり、== のところが、= だったり、Error は様々。
File not found のときは、データを読み込むコードを再実行。それでもうまくいかないときは、data フォルダに、データ（***.csv）が入っているかを確認、なければ、data フォルダがあることを確認、なければ作成して読み込みを実行。
Object not found のエラーがでたときには、上から順に、Run （Code Chunk の右上の三角印を押して実行）。または、エラー箇所にカーソルを置き、Run ボタンの右の三角から、Run All Chunks Above を選択して実行。
実行できていても、結果が見えないことがある。そのときは、Code Chunk の下にある、山二つの記号を押す。結果の表示、非表示の切り替え。
Error message を読む。エラーが出た Code Chunk と、Error message を、ChatGPT や、Google Bard, Google Search に入れ解決方法を探る。
Posit.Cloud で、Error 137 が出た時はメモリーオーバーなので、例など不要な部分を削除。ただし、一番上の、YAML は、かならず残すこと。

参考文献

課題1：解説リンク・解答例リンク
探索的データ分析（課題2の解説にかえて）in progress リンク

「みんなのデータサイエンス - Data Science for All」[はじめてのデータサイエンス]
Posit Recipes（旧 Posit Primers）: The Basics 対話型の演習サイトの最初 [Link]
Posit Cheat Sheet. 早見表です。印刷して使うために、PDF も提供しています。[Site Link]
DataCamp Cheat Sheet: Tidyverse for Biginners. データサイエンスの教育をしている会社の早見表の一つです。基本が簡単にまとまっています。[Link]
箱ひげ図の見方 Link、外れ値検出のある箱ひげ図 Link
Video: Histogram, Boxplot

準備