第6週 探索的データ分析 - EDA, Part I

演習

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

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

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

データ情報

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

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

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

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

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

確認作業

• 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", "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_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 = "日本の一人当たりの二酸化炭素排出量")

気づいたこと・疑問

• 9トン前後で推移している。
• 2020年に落ちているのは、コロナのせいかもしれないので、このあとの変化が重要。

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

気づいたこと・疑問

• 0.25% 程度の幅で推移している。
• 日本の国土が、380,000 平方キロとすると、0.25% は、950 平方キロなので、だいたい、30km x 30km の森林が増減しているということか。やはり、原因を知りたい。

b. 選択した国・地域

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

気づいたこと・疑問

• 中国が急激に増加している。1990年の三倍程度。
• ロシアは、1993年のソ連崩壊以後、しばらく減少。

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

気づいたこと・疑問

• 中国とインドは、森林面積も増加している。
• ブラジルの1990年から2010年の間の、森林面積の減少は10%程度あるようだ。

二つのデータの散布図

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

a. すべての値の散布図

df_w6 |> ggplot(aes(forest, co2pcap, col = region)) + 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

実習

例：二酸化炭素排出量と[選択した指標]

データ

データ情報

• データ１：一人当たりの二酸化炭素排出量 (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.

• データ２：名前、コード、変数名、リンク

• 概要：

データの取得

準備

library(tidyverse)
library(WDI)

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

データの名前は、変えたほうがよいので、例でも、df_w6eda_1 や、df_w6_1 に変えてあります。

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

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

write_csv(df_w6eda_1, "data/w6eda_1.csv")

df_w6eda_1 <- read_csv("data/w6eda_1.csv")

データの確認

df_w6eda_1

df_w6eda_1

Error: object 'df_w6eda_1' not found

str(df_w6eda_1)

変数の選択（selecting）

df_w6_1 <- df_w6eda_1 |> 
  select(country, iso2c, year, co2pcap, forest, region, income)
df_w6_1

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

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

国と地域

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_w6eda_1 |> filter(iso2c %in% REGION) |> distinct(country, iso2c)

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

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

分析する国のリスト

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

経年変化

a. 日本

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

気づいたこと・疑問

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

気づいたこと・疑問

b. 選択した国・地域

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

気づいたこと・疑問

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

気づいたこと・疑問

二つのデータの散布図

必要に応じて log10 スケール

a. すべての値の散布図

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

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

df_w6_1 |> 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_1 |> 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_1 |> 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)

気づいたこと・疑問

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

気づいたこと・疑問