GES 001 演習6

2024年1月31日

準備

（自分の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追加提出ボックスをMoodle内に作成しました。そこに、課題1の解説や、11個の指標についての解答例もつけました。確認してください。

第6週

01/25(TH)　気候変動問題の原因：気候変動と経済活動1

　　　　　　気候変動問題の原因：気候変動と経済活2

第6週、第7週の講義では、気候変動が我々の生活にもたらす影響と対策を議論します。

基本的にIPCCの第6次評価報告書（https://www.data.jma.go.jp/cpdinfo/ipcc/index.html ）

をテキストに使っています。

01/30(TU)　Rでデータサイエンス6：気候変動　 [Main]・[授業]

Environment-Climate　

CO2 emissions (metric tons per capita) ：EN.ATM.CO2E.PC [Link]

Forest area (% of land area)：AG.LND.FRST.ZS [Link]

Renewable electricity output (% of total electricity output)：EG.ELC.RNEW.ZS [Link]

Electricity production from oil, gas and coal sources (% of total)：EG.ELC.FOSL.ZS [Link]

Electricity production from nuclear sources (% of total)：EG.ELC.NUCL.ZS [Link]

IEA:　Energy Statistics Data Browser　データ元

https://www.iea.org/data-and-statistics/data-tools/energy-statistics-data-browser?country=WORLD&fuel=Energy%20supply&indicator=ElecGenByFuel

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

表題 Title
概要 Abstract
データ情報：データ名、データコード、変数名、データ概要
データの取得 - ダウンロードして保存しておく
- 準備：library(tidyverse); library(WDI)
- WDI(indicator(ed_exp = "SE.XPD.TOTL.GD.ZS")) 一つの指標
- データの確認：データが十分あるか、国名リスト、地域名リスト
分析する国のリスト：南部アフリカ関税同盟、国や地域を選択

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

分析・視覚化
- 各年毎のデータの数の棒グラフ
- 経年変化
  - 日本（身近な馴染みのあるデータを確認）
  - 南部アフリカ関税同盟、選択した国・地域
- 分布
  - ヒストグラム（年を指定した度数分布）
  - 値が大きい10カ国とその値の棒グラフ
  - 値が小さい10カ国とその値の棒グラフ
各項目に、気づいたこと、疑問などを記録

演習 1月31日（火）

内容

データ情報

CO2 emissions (metric tons per capita) ：EN.ATM.CO2E.PC [Link]
Forest area (% of land area)：AG.LND.FRST.ZS [Link]
Renewable electricity output (% of total electricity output)：EG.ELC.RNEW.ZS [Link]
Electricity production from oil, gas and coal sources (% of total)：EG.ELC.FOSL.ZS [Link]
Electricity production from nuclear sources (% of total)：EG.ELC.NUCL.ZS [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_environment <- WDI(
  indicator = c(co2pcap = "EN.ATM.CO2E.PC",
                forest = "AG.LND.FRST.ZS",
                renewable = "EG.ELC.RNEW.ZS",
                fossil = "EG.ELC.FOSL.ZS",
                nuclear = "EG.ELC.NUCL.ZS"
                ), extra = TRUE)

２回目からは、data から読み込めるようにしておく ファイル (Rmd) の保存場所に data フォルダがあることを確認

write_csv(df_environment, "data/environment.csv")

df_environment <- read_csv("data/environment.csv")

## Rows: 16758 Columns: 17
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (7): country, iso2c, iso3c, region, capital, income, lending
## dbl  (8): year, co2pcap, forest, renewable, fossil, nuclear, longitude, lati...
## 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_environment

地域の `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_environment |> filter(iso2c %in% REGION) |> distinct(country, iso2c)

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

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

変数の選択（selecting）

df_env <- df_environment |> 
  select(country, iso2c, year, co2pcap, forest, 
         renewable, fossil, nuclear, region, income)
df_env

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

df_env |> filter(!(iso2c %in% REGION)) |> # 地域以外（国のみ選択）
  drop_na(co2pcap, forest, renewable, fossil, nuclear) |> 
  select(co2pcap, forest, renewable, fossil, nuclear) |> cor()

##              co2pcap     forest  renewable     fossil    nuclear
## co2pcap    1.0000000 -0.1832374 -0.4138267  0.3220169  0.1369711
## forest    -0.1832374  1.0000000  0.3683177 -0.4234481  0.1131410
## renewable -0.4138267  0.3683177  1.0000000 -0.8784240 -0.1851655
## fossil     0.3220169 -0.4234481 -0.8784240  1.0000000 -0.2354258
## nuclear    0.1369711  0.1131410 -0.1851655 -0.2354258  1.0000000

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

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

ファイルリンク

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

探索的データ分析（EDA） -
- 練習と一つ目の実習課題（w6eda.Rmd） [リンク], [Rmd]
- 二つ目以降の実習課題：w6eda1.Rmd [リンク], [Rmd]

R Notebook のソースファイル（Rmd）を取得方法

Rmd のリンクをクリックし Raw の横の Copy a raw file からコピー、新規 RMarkdown ファイルを PositCloud または RStudio のプロジェクト内に作成しペースト

演習

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

各年毎のデータの数の棒グラフ
国または地域を選択
それぞれの経年変化を表す折れ線グラフ
1. 日本
2. 選択した国または地域
二つのデータの散布図-必要に応じて log10 スケールを用いる
1. すべての値の散布図
2. NA ではない値の散布図、近似（回帰）直線を表示
3. 地域を除き国のみの散布図、近似（回帰）直線を表示
4. 最近の年を選択し、地域を除き国のもの散布図、近似（回帰）直線を表示

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

確認作業

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

途中でのエラー

入力したときには、例を参照して、スペルなどを確認してください。全角になっていると問題がおきます。() がペアでマッチしているか、確認してください。
引用符が入っていなかったり、== のところが、= だったり、いろいろな可能性があります。Error message を読むこともたいせつです。エラーがでた、Code Chunk と、Error message を、ChatGPT や、Google Bard, Google Search に入れると、解決方法を教えてくれることもあります。
File not found の、エラーがでたときには、上から順に、Run （Code Chunk の右上の三角印を押して実行）してみてください。または、エラーが出たところに、カーソルを置き、上の、Run ボタンの右の三角から、Run All Chunks Above を選択すると、そこまでのすべての　Code Chunk を実行してくれます。
上の方法でうまくいかないときは、data フォルダに、データ（***.csv）が入っているかを確認、なければ、data フォルダがあることを確認して、最初のデータ読み込みのところを実行してみてください。
実行できていても、結果が見えないこともあります。そのときは、Code Chunk の下にある、山二つの記号を押してみてください。これは、結果を表示、非表示にします。それが原因で隠れている場合があります。

データ

データ情報

データ１：一人当たりの二酸化炭素排出量 (CO2 emissions (metric tons per capita))、“EN.ATM.CO2E.PC”、co2pcap [Link]
概要：Carbon dioxide emissions are those stemming from the burning of fossil fuels and the manufacture of cement. They include carbon dioxide produced during consumption of solid, liquid, and gas fuels and gas flaring.
データ２：森林面積（％）(Forest area (% of land area))、“AG.LND.FRST.ZS”、forest [Link]
概要：Forest area is land under natural or planted stands of trees of at least 5 meters in situ, whether productive or not, and excludes tree stands in agricultural production systems (for example, in fruit plantations and agroforestry systems) and trees in urban parks and gardens.

データの取得

準備

library(tidyverse)
library(WDI)

データのダウンロードと保存：コードと変数名を指定。

df_w6eda <- WDI(indicator = c(co2pcap = "EN.ATM.CO2E.PC",
                              forest = "AG.LND.FRST.ZS"),
                extra = TRUE)

２回目からは、data から読み込めるようにしておくファイル (Rmd) の保存場所に data フォルダがあることを確認

write_csv(df_w6eda, "data/w6eda.csv")

df_w6eda <- read_csv("data/w6eda.csv")

## Rows: 16758 Columns: 14
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (7): country, iso2c, iso3c, region, capital, income, lending
## dbl  (5): year, co2pcap, forest, longitude, latitude
## 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_w6eda

str(df_w6eda)

## spc_tbl_ [16,758 × 14] (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" ...
##  $ co2pcap    : num [1:16758] 0.2837 NA 0.0898 0.2981 0.0888 ...
##  $ forest     : num [1:16758] 1.85 NA 1.85 1.85 1.85 ...
##  $ 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 = ""),
##   ..   co2pcap = col_double(),
##   ..   forest = 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_w6 <- df_w6eda |> 
  select(country, iso2c, year, co2pcap, forest, region, income)
df_w6

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

df_w6eda |> drop_na(co2pcap, forest) |>
  ggplot(aes(year)) + geom_bar()

国と地域

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

REGION <- c("1A", "1W", "4E", "7E", "8S", "B8", "EU", "F1", "OE", "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", "ZF", 
"ZG", "ZH", "ZI", "ZJ", "ZQ", "ZT")

地域のリストを表示

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

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

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

分析する国のリスト

BRICS を選択します。

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

経年変化

a. 日本

df_w6 |> drop_na(co2pcap) |> filter(country == "Japan") |>
  ggplot(aes(year, co2pcap)) + geom_line() + 
  labs(title = "日本の一人当たりの二酸化炭素排出量")

df_w6 |> drop_na(forest) |> filter(country == "Japan") |>
  ggplot(aes(year, forest)) + geom_line() + 
  labs(title = "日本の森林面積（％）")

b. 選択した国・地域

df_w6 |> drop_na(co2pcap) |> filter(country %in% BRICS) |>
  ggplot(aes(year, co2pcap, col = country)) + geom_line() + 
  labs(title = "BRICS の一人当たりの二酸化炭素排出量")

df_w6 |> drop_na(forest) |> filter(country %in% BRICS) |>
  ggplot(aes(year, forest, col = country)) + geom_line() + 
  labs(title = "BRICSの森林面積（％）")

二つのデータの散布図

必要に応じて log10 スケール (+ scale_y_log10)

a. すべての値の散布図

df_w6 |> ggplot(aes(forest, co2pcap, col = region)) + geom_point()

## Warning: Removed 9494 rows containing missing values (`geom_point()`).

df_w6 |> filter(co2pcap >0) |> ggplot(aes(forest, co2pcap, col = region)) + 
  geom_point() + scale_y_log10()

## Warning: Removed 139 rows containing missing values (`geom_point()`).

b. NA ではない値の散布図、近似（回帰）直線を表示

df_w6 |> drop_na(co2pcap, forest) |> 
  ggplot(aes(forest, co2pcap)) + geom_point(aes(col = region)) +
  geom_smooth(formula = 'y~x', method = "lm", se = FALSE)

c. 地域を除き国のみの散布図、近似（回帰）直線を表示

df_w6 |> filter(!(iso2c %in% REGION)) |> drop_na(co2pcap, forest) |> 
  ggplot(aes(forest, co2pcap)) + geom_point(aes(col = region)) +
  geom_smooth(formula = 'y~x', method = "lm", se = FALSE)

d. 最近の年を選択し、地域を除き国のもの散布図、近似（回帰）直線を表示

df_w6 |> filter(!(iso2c %in% REGION)) |> filter(year == 2020) |> drop_na(co2pcap, forest) |> 
  ggplot(aes(forest, co2pcap)) + geom_point(aes(col = region)) +
  geom_smooth(formula = 'y~x', method = "lm", se = FALSE)

気づいたこと・疑問

森林面積が10％よりも少ない国が多い
弱い負の相関があるようだ
授業で求めた相関係数は、-0.1832374。これは、四つの指標すべての値があり（NA ではなく）国地域も分けず、すべての年のデータでの相関をとったもの。傾向はある程度わかる。

df_w6 |> filter(!(iso2c %in% REGION)) |> filter(year == 2020) |> drop_na(co2pcap, forest) |> 
  select(co2pcap, forest) |> cor()

##             co2pcap      forest
## co2pcap  1.00000000 -0.09914706
## forest  -0.09914706  1.00000000

参考

df_w6 |> filter(!(iso2c %in% REGION)) |> filter(year == 2020) |> drop_na(co2pcap, forest) |> 
  select(co2pcap, forest) |> lm(co2pcap ~ forest, data = _) |> summary()

Residuals: 直線よりどのぐらい大き値や小さい値があるか
Coefficients: Estimate の　(Intercept) は直線の y-切片、forest は、傾き
Residual Standard Error: 直線からこの値の幅の間に、おおよそ 2/3 程度（68.2%）の点が入っている
Multiple R-squared: 直線がどの程度近似しているか、0 と 1 の間で、相関係数の二乗
p-value: かなり小さいと（たとえば 0.01 未満）相関が正か負かがほぼ確実といえる

## 
## Call:
## lm(formula = co2pcap ~ forest, data = select(drop_na(filter(filter(df_w6, 
##     !(iso2c %in% REGION)), year == 2020), co2pcap, forest), co2pcap, 
##     forest))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.306 -2.958 -1.105  1.216 27.315 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.41204    0.56135    7.86 2.84e-13 ***
## forest      -0.01897    0.01385   -1.37    0.172    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.598 on 189 degrees of freedom
## Multiple R-squared:  0.00983,    Adjusted R-squared:  0.004591 
## F-statistic: 1.876 on 1 and 189 DF,  p-value: 0.1724

library(plotly)
df_w6_chart_f <- df_w6 |> filter(!(iso2c %in% REGION), year == 2020) |> drop_na(co2pcap, forest) |> 
  ggplot(aes(forest, co2pcap, text = country, colour = region)) + geom_point() +
  labs(title = "Forest (%) vs CO2 per Capita (tons)") + theme(legend.position = "none")
df_w6_chart_f |> ggplotly()

演習の内容と課題

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

探索的データ分析（EDA） -
- 練習と一つ目の課題（w6eda.Rmd） [リンク], [Rmd]
- 二つ目以降の課題：w6eda1.Rmd [リンク], [Rmd]

参考文献

課題1：解説リンク
課題1：解答例リンク

「みんなのデータサイエンス - Data Science for All」[はじめてのデータサイエンス]
- 導入として、GDP（国内総生産）のデータを使って説明しています。
Posit Recipes（旧 Posit Primers）: The Basics 対話型の演習サイトの最初 [Link]
Posit Cheat Sheet. 早見表です。印刷して使うために、PDF も提供しています。[Site Link]
DataCamp Cheat Sheet: Tidyverse for Biginners. データサイエンスの教育をしている会社の早見表の一つです。基本が簡単にまとまっています。[Link]

準備

ファイル

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

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

課題1について

第6週

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

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

演習 1月31日（火）

内容

データ情報

データの取得

準備

データのダウンロードと保存

データの確認

地域の iso2c コード

変数の選択（selecting）

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

ファイルリンク

演習

確認作業

途中でのエラー

データ

データ情報

データの取得

準備

データの確認

変数の選択（selecting）

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

国と地域

地域のリストを表示

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

分析する国のリスト

経年変化

a. 日本

b. 選択した国・地域

二つのデータの散布図

a. すべての値の散布図

b. NA ではない値の散布図、近似（回帰）直線を表示

c. 地域を除き国のみの散布図、近似（回帰）直線を表示

d. 最近の年を選択し、地域を除き国のもの散布図、近似（回帰）直線を表示

参考

演習の内容と課題

参考文献

地域の `iso2c` コード