探索的データ分析（Exploratory Data Analysis）

1 ワークフロー

• 表題 Title
• 概要 Abstract
• データ情報：データ名、データコード、変数名、データ概要
• データの取得 - ダウンロードして保存しておく
  • 準備：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 収入レベル
  • データの確認：データが十分あるか、国名リスト、地域名リスト
• 分析する国のリスト：南部アフリカ関税同盟、国や地域を選択
• 分析・視覚化
  • 各年毎のデータの数の棒グラフ
  • 経年変化
    • 日本（身近な馴染みのあるデータを確認）
    • 南部アフリカ関税同盟、選択した国・地域
  • 分布
    • ヒストグラム（年を指定した度数分布）
    • 値が大きい10カ国とその値の棒グラフ
    • 値が小さい10カ国とその値の棒グラフ
  • 相関
    • 散布図（scatter plot）x-軸と、y-軸に、それぞれの変数をとる
    • 相関係数
• 各項目に、気づいたこと、疑問などを記録

2 準備

Posit.Cloud または RStudio で R を使うことを想定して説明と注意点を書きます。

2.1 R Notebook の作成

再現性と記録、さらには、コミュニケーションのために、R Notebook を作成します。R Notebook は、R Markdown の一形式です。

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

1. Posit Cloud のときは、まず、Login し、プロジェクトに移動（クラスでは、共有プロジェクト intro2rj を使用）のプロジェクトに入り、 ファイルから、ges001 のフォルダーを選択して、移動します。この中に、data フォルダが作成されていることを確認し、このges001 に、新しく作成したファイルを保存します。新しく作成したファイルの入っているフォルダーの中に、data フォルダがあることが大切です。そこに、データを書き込みます。

2. RStudio の場合には、まずは、新しい Project を作成します。File > New Project から作成します。すでに、作成してある場合は、それを、Open Project や、Recent Project から開きます。その中に、新しいファイルを作成します。作成したフォルダーに、data フォルダがあることを確認してください。新しく作成したファイルの入っているフォルダーの中に、data フォルダがあることが大切です。そこに、データを書き込みます。

• 　Rmd ファイルを作成したら、すぐに、ID と名前を書きます。タイトルも決まっていたら書きます。

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

2.1.2.1 確認作業

• Run All で、最後に、もう一度、すべてを、実行。

• エラーが修正できない場合には、その次のコードチャンクから、Run Below で実行。

• Preview の右のギアマークから、Preview in Widow を選択し、View in Browser で、いつも、使っている、Web Browser で、最後まで、グラフなどが出ているか確認。

2.1.2.2 Posit.Cloud のファイルダウンロード

RStudio の場合には、自分の PC に作成したファイルがありますから、問題ないと思いますが、Posit Cloud で作成した場合には、自分の PC にダウンロードする必要があります。

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

2.2 必要なパッケージの読み込み

このコースでは、この二つのパッケージ（R の機能を拡張するもの）を読み込みます。初めて使う時は、Tools > Install Packages からインストールします。

library(tidyverse)
library(WDI)

3 データ

• 世界開発指標（World Development Indicators）Link からデータを選択します。

• 授業で使ったデータのリスト：リンク

3.1 データ情報

データ名、データID（WDI コード）、データ概要を取得して書きます。

リンクを残しておくとあとで便利です。サイトの Details に概要も、データ ID もあります。日本語にしておくのも良いでしょう。

このときに、短い変数名も決めておきます。

3.1.1 例

データ１：一人当たりの二酸化炭素排出量（トン）CO2 emissions (metric tons per capita) ：EN.ATM.CO2E.PC [Link]

データ１変数名：co2

データ１概要：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.

データ２：化石燃料での電気生成の割合（%）Electricity production from oil, gas and coal sources (% of total)：EG.ELC.FOSL.ZS [Link]

データ２変数名：fossil

データ２概要：Sources of electricity refer to the inputs used to generate electricity. Oil refers to crude oil and petroleum products. Gas refers to natural gas but excludes natural gas liquids. Coal refers to all coal and brown coal, both primary (including hard coal and lignite-brown coal) and derived fuels (including patent fuel, coke oven coke, gas coke, coke oven gas, and blast furnace gas). Peat is also included in this category.

3.2 データの読み込み（Importing）

WDI で読み込みます。

• 読み込むときには、指標 ID と、変数名が必要です。
• 読み込んだデータに名前をつけます。

3.2.1 例

この例では、二つの指標を、読み込みます。

df_co2_fossil <- WDI(indicator = c(co2 = "EN.ATM.CO2E.PC",
                     fossil = "EG.ELC.FOSL.ZS"), extra = TRUE)
df_co2_fossil

• 一つ目のデータの指標（ID）は、EN.ATM.CO2E.PC、変数名は、co2 にしました。二つ目のデータの指標（ID）は、EG.ELC.FOSL.ZS、変数名は、fossil にしました。

• 指標 (ID) は WDI のサイトから得たものを使わなければいけませんが、変数名は自分で決められます。 できるだけ、分かりやすい名前にしておくのが良いでしょう。データ名と、変数名を区別する為、わたしは、データは、df_ で始まる名前にしています。

• 最後にextra = TRUE としておくと、国ごとの地域、所得レベルなどの情報を一緒に読み込んでくれます。必要がなければ、次のようになります。

df_co2_fossil_short <- WDI(indicator = c(co2 = "EN.ATM.CO2E.PC",
                     fossil = "EG.ELC.FOSL.ZS"))
df_co2_fossil_short

-「機械は思った通りに動かない。書いた通りに動く。」と言われます。引用符や、かっこなど、マッチするように、書いてください。

3.2.2 WDI でのデータの読み込みのコードについて

実は、次の命令でも、データをとってきます。

WDI()

Help をみると次のように書かれています。

WDI(
  country = "all",
  indicator = "NY.GDP.PCAP.KD",
  start = 1960,
  end = NULL,
  extra = FALSE,
  cache = NULL,
  latest = NULL,
  language = "en"
)

これらの変数を指定できますが、右にあるのが、初期値です。一つ一つの説明も、Help ページに書いてあります。

WDI() だと、GDP (PPP) のすべての国のデータとダウンロードします。使う為には、これに名前をつける必要があるので、上の例では、df_co2_fossil とか、df_co2_fossil_short としてあります。 最後の行に、これらを加えておくことで、内容を表示することも同時にしてくれます。

3.3 データの保存

ダウンロードしたデータは、data フォルダに保存しておくことをお勧めします。世界銀行は、ときどき、データをアップデータするために、サービスが停止することがあります。また、保存したものを読み込むことで、毎回ダウンロードしないで済むメリットもあります。

3.3.1 例

読み込んだデータ df_co2_fossil を、CSV（カンマ区切り形式）で、data フォルダーの中に、co2_fossil.csv という名前で保存します。このときに、data フォールだがないと、保存できません。R Notebook があるフォルダーの中に、data フォルダがあることが必要です。data フォルダがないというエラーメッセージが出たら、最初に説明した、ファイルを作る時の部分を読み返してください。間違っていたら、ここで、Save as で適切なところに、保存しなおしてください。

write_csv(df_co2_fossil, "data/co2_fossil.csv")

df_co2_fossil <- read_csv("data/co2_fossil.csv")

Rows: 16758 Columns: 14── Column specification ────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (7): country, iso2c, iso3c, region, capital, income, lending
dbl  (5): year, co2, fossil, 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.

write_csv(df_co2_fossil_short, "data/co2_fossil_short.csv")

df_co2_fossil_short <- read_csv("data/co2_fossil_short.csv")

Rows: 16758 Columns: 6── Column specification ────────────────────────────────────────────────────────────────
Delimiter: ","
chr (3): country, iso2c, iso3c
dbl (3): year, co2, fossil
ℹ 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.

3.4 データの確認

データを確認します。どのような列（変数）があるかを確認しますが、特に、最初に定めた、指標の変数が、ある列にはっていることを確認することも大切です。特に、extra = TRUE として読み込むと、さまざまな変数を読み込んでいますから、確認してください。

3.4.1 例

データ全体の確認。右上の三角印から、隠れている変数をみることができ、下の NEXT や、番号から、次のページへと移動し、10データずつ確認できるようになっています。

head(df_co2_fossil) とすると、最初の６行だけ、tail(df_co2_fossil) とすると、最後の、6行を確認することができます。

df_co2_fossil

データの変数（列）の確認には、str() も便利です。下のコードは、df_co2_fossil |> str() としても同じです。|> はパイプオペレータと呼ばれます。かっこの中の第一引数として、一つ前の出力を代入します。

str(df_co2_fossil)

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-19" "2023-12-19" ...
 $ co2        : num [1:16758] 0.2837 NA 0.0898 0.2981 0.0888 ...
 $ fossil     : 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 = ""),
  ..   co2 = col_double(),
  ..   fossil = col_double(),
  ..   region = col_character(),
  ..   capital = col_character(),
  ..   longitude = col_double(),
  ..   latitude = col_double(),
  ..   income = col_character(),
  ..   lending = col_character()
  .. )
 - attr(*, "problems")=<externalptr> 

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

世界銀行では、1960年から昨年までのデータを作成していますが、すべてあるわけではありません。また、国や地域によっては、データがない場合があります。最新のデータは何年のものかも含めて、調べておく必要があります。

3.5.1 例

この棒グラフでは、1990年から2015年までのデータが180程度、すなわち、180の国の地域については、両方の知っ方についてのデータがあることがわかります。

df_co2_fossil |> drop_na(co2, fossil) |> 
  ggplot(aes(year)) + geom_bar() + labs(title = "各年毎のデータの数：co2, fossil ともに存在する国・地域の数")

各、変数ごとに調べたい時は、long table を使います。少し高級ですね。

この例では、fossil のデータは、1960 年からありますが、co2 は、1990 年から。また、co2 は、2020 年までのデータがありますが、fossil は、2015 年までのデータしかないことがわかります。

pivot_longer(c(co2, fossil)) では、これら二つの指標を、name という変数に、それぞれの値は、value に入れます。縦長のデータということで、説明しました。（演習7）

df_co2_fossil |> pivot_longer(c(co2, fossil)) |> drop_na(value) |> 
  ggplot(aes(year, fill = name)) + geom_bar()+ 
  labs(title = "各年毎のデータの数：co2, fossil それぞれが存在する国・地域の数")

4 変形（transforming）

データの必要な部分だけを取り出して、分析します。

4.1 基本的な変形

• 行の選択：filter(country %in% c("Japan", "China")), distinct(country), drop_na(value)

• 列の選択：select(country, iso2c, year, var1, var2, region)

4.2 縦長データ・横長データ

• 横並びの変数を縦並びに：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(c(electricity_urban, electricity_rural)): cols = は省略可能

  • pivot_longer(cols = electricity_urban:sanitation_rural))：electricity_urban の列から、sanitation_rural の列までのすべての変数を指定

• 活用：pivot_longer(cols = c(var1, var2))|> ggplot(aes(year, value, col = name) + geom_line()

• 縦並びの変数を横並びに：列名を指定する変数と、値の変数を指定し、横並びにします。

  • `pivot_wider(names_from = var1, values_from = var2)

  • `pivot_wider(names_from = name, values_from = value): pivot_longer の逆です。

4.3 変数の選択（selecting）

extra = TRUE として、データを入手すると、使わない変数もたくさんあるので、見やすいように、必要な変数だけにします。必要なものを、落とさないように、気をつけます。

4.3.1 例

わたしは、少し短めの名前を使っています。結果を、すぐ表示するようにして、そこで、必要な変数が、のこていることを確認します。

df_co2fos <- df_co2_fossil |> 
  select(country, iso2c, year, co2, fossil, region, income)
df_co2fos

4.4 縦長の表（Long Table）に変換

上の説明を、確認してください。たいせつなのは、これをどのように活用するかですから、このあとに、df_co2fos_long を使って、図を描いているものを、確認して、df_co2fos との違いを比較してください。このあとも、df_co2fos を使って描画をしている場合もあります。どちらでも、使えるものと、一方を使うと、他方を使うのと比較して、非常に、単純になる場合があることを確認していってください。

4.4.1 例

co2 と fossil の列（変数）を、変数を name に、値を value に入れて、縦に並べます。

df_co2fos_long <- df_co2fos |> 
  pivot_longer(cols = c(co2,fossil))
df_co2fos_long

4.5 国と地域

主として分析する国や地域のリストを作成します。最初に書いてありますが、このあとの分析で、国や、地域を加えたり、変更することもあるでしょう。全体の構成としては、視覚化などを始める前に、リストを作成しておくことが必要です。定義がされていないと、エラーが生じるからです。最後に、Run All などをすればわかりますが、上から順に実行するので、たとえば、ASEAN などの分析をするときに、それが、最初の方で定義されていなければ、エラーとなります。途中では、その都度定義すれば、それを、理解して実行しますが、完成形にするときには、実際に使うコードより前に定義しておかないといけません。

国や地域を選択するには、国名 country または、iso2c などの国のコードが必要です。正確でないといけませんし、iso2c などをすべて覚えているわけではないので、それらを表示させ、そこから、コピー・ペーストするのが安全です。

国名は一般的に長いので、iso2c も便利です。指定する時は、たとえば、asean に国名が、ASEAN に iso2c コードが入っていれば、

• filter(country %in$ asean)

• filter(iso2c %in% ASEAN)

選択するときの、条件が変わりますから、注意してください。

4.5.1 地域の iso2c コード

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

• 地域のみ filter(country %in% REGION)
• 国のみ filter(!(country %in% REGION))

! は否定です。

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")

4.5.1.1 地域のリストを表示

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

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

df_co2_fossil とそれから変数を選択した、df_co2fos は、extra=TRUE として、追加情報がありますから、地域情報（region）と 所得レベル（income）も加えておきます。

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

df_co2_fossil_short では、extra=TRUE は指定せず（初期設定通り extra = FALSE）としてありますから、上の、df_co2fos を、df_co2_fossil_short に置き換えるとエラーが出ます。region や、income の列がないからです。その場合は、次のようにします。

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

存在しない変数についてなにかをしようとすると Error になります。そのときは、そこで使っているデータに実際に、その変数があるかどうか確認してください。この場合なら、以下のようにすれば、region や、income という変数がないことは確認できますね。

df_co2_fossil_short

4.5.2 国や地域の選択

授業で紹介したものを次に書いておきます。

asean <- c("Brunei Darussalam", "Cambodia", "Indonesia", "Lao PDR", "Malaysia", 
"Myanmar", "Philippines", "Singapore", "Viet Nam")

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")

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

4.5.3 地域（region）と、所得レベル

extra = TRUE として、データを入手すると、地域名と、所得レベルも一緒に得られます。それらは、どのようなものか確認してみます。distinct(region) は、異なるもののみ取得、pull() は、表形式ではなく、ベクトル（データのならび）で取得する為に使っています。

df_co2fos |> distinct(region) |> pull()

[1] "South Asia"                 "Aggregates"                
[3] "Europe & Central Asia"      "Middle East & North Africa"
[5] "East Asia & Pacific"        "Sub-Saharan Africa"        
[7] "Latin America & Caribbean"  "North America"             
[9] NA                          

これも、定義しておきます。REGION はすでに使ってありますから、混乱をさけるため、小文字にしました。

regions <- c("South Asia", "Europe & Central Asia", "Middle East & North Africa", 
"East Asia & Pacific", "Sub-Saharan Africa", "Latin America & Caribbean", 
"North America")

df_co2fos |> distinct(income) |> pull()

[1] "Low income"          "Aggregates"          "Upper middle income"
[4] "Lower middle income" "High income"         NA                   
[7] "Not classified"     

こちらは、ある程度順番も重要ですから、次のようにしておきます。

INCOME <- c("High income", "Upper middle income", "Lower middle income", 
"Low income")

4.5.4 分析する国や地域のリスト

あなたが分析したい国や地域のリストをまずは書いておくと、思いつきで変更するのではなく、首尾一貫した分析ができます。

最初の確認のためには、わたしは、Japan (iso2c では、JP) か、World (iso2c では、1W) をつかうことにしています。身近な、または、全体から始める為です。分析を始めたら、CHOSEN_COUNTRIES は変更した方が良い場合もありますし、CHOSEN_COUNTRIES_1, CHOSEN_COUNTRIES_2 などとして、幾つものグループについて分析するのもよいでしょう。もっと短い名前にしてもよいですが、途中まで、入力すればあとは自動的に表示されますから、矢印と、Tab を使えば、簡単に入力できます。

CHOSEN_COUNTRY <- ""
CHOSEN_COUNTRIES <- c()

5 分析・視覚化

一般的な形で書く時は、今後、データは、df で表します。data frame の意味です。R では、基本的なデータ形式です。df_co2fos などとした、df もそこからとったものです。

視覚化は、探索的データ分析の核とも言えるものです。視覚化して、それから、見えるものを通して、考えます。気づいたこと、疑問に思ったことを、書き留めておく。そして、次の視覚化で、さらに、理解を深めていくといった感じです。また、ちがった、データをみることもあるかと思いますが、ここでは、限られたデータをみる、基本的な方法を見ていきます。

基本は、二つ。

1. 一つのデータの、分布をみる。
2. 二つのデータの相関をみる。

データは、数値データ（numerical data）ばかりでなく、文字からなるデータ、通常は、カテゴリー・データ（categorical data、日本語では、名目データということもあります）の場合もあります。

5.1 さまざまな視覚化（Visualizing）

例として、中心的に使うデータは、二つです。

• df_co2fos: co2 (Co2 per Capita) の値と、化石燃料による発電の割合を表す fossil が変数に含まれ、year 以外に、region や、income のデータも入っています。

• df_co2fos_long: 縦長の表で、name に、co2 と fossil を含み、value に その値が入っています。データをよく確認してください。上の方にあります。

基本的には、データを絞り込み、そのあとで、ggplot2 を使います。以下のようなフォーマットになります。

df_co2fos |> filter などで必要な行を絞り込む　|> drop_na(変数) などで、NA を削除 |>
  ggplot(aes(x, y)) + geom_***()

aes() の部分を mapping と言います。実は、mapping = aes(x, y) と書くのが正式で、その　mapping を省略しています。詳細は、geom_***() のいろいろな描写方法によって異なりますので、それぞれの、Help（右下の窓枠（pane）の Help タブの検索窓に入れれば、説明や、パラメタ（オプションとして書き込む内容）や、その初期値などの情報が得られます）を確認してください。

データの最初の6行を表示させるのは、head() でしたが、head と Help に入れてみてください。かなり詳しく書かれていますから、すべて理解するのは難しいかもしれませんね。視覚化の最初は、折れ線グラフから始めますが、それには、geom_line を使います。これもみてみてください。

5.2 経年変化

year を x 軸に、変数の値を　y 軸に取ると、経年変化がわかります。日本や、世界でみてみます。drop_na(var) を使うことで、空白部分がなくなり、Warning も減りますが、違う指標を比較する時は、これを入れないのも手ですし、同じ基準にするために、drop_na(var1, var2) としておく手もあります。

基本は、年毎に、数値変数 var の点を、結びます。

df |> ggplot(aes(year, var)) + geom_line()

次のように書くことも可能です。

df |> ggplot() +geom_line(aes(year, var))

df |> ggplot(aes(year)) + geom_line(var1) + geom_line(var2)

なども可能です。

例をみていきましょう。一つの国・地域

df_co2_fos |> filter(country == "国名") |> drop_na(変数) |>
  ggplot(aes(year, 変数)) + geom_line()

いくつかの国・地域

df_co2_fos |> filter(country %in% c("国名1", "国名2", ... , "国名n")) |> drop_na(変数) |>
  ggplot(aes(year, 変数, col = country)) + geom_line()

iso2c で、国・地域を選択することも可能です。

df_co2_fos |> filter(iso2c %in% c("二文字コード1", "二文字コード2", ... , "二文字コードn")) |>
  drop_na(変数) |>
  ggplot(aes(year, 変数, col = country)) + geom_line()

最初分析中は、省略してもよいですが、なるべく、

• labs(title = “グラフのタイトル”)

として、タイトル名を入れましょう。x 軸、y 軸のラベルを変更する時は、

• labs(title = “グラフのタイトル”, x = “年”, y = “変数の名前”)

とします。

改行するときには、|> の後、+、(、, や、= の後で改行します。なんとなく、続くことが明らかなところで改行するということですね。

5.2.1 例

5.2.1.1 国または地域の一つの変数に関する経年変化

一番、簡単な形は、

df_co2fos |> filter(country == "Japan") |> ggplot(aes(year, co2)) + geom_line() +
  labs(title = "日本の一人当たりの二酸化炭素排出量：NA も含む")

ただ、これだと、最初の方はデータがないので、グラフが途中から始まります。そこで、drop_na(co2) を入れて、co2 の値が、NA ではないところだけ表示します。

df_co2fos |> filter(country == "Japan") |> drop_na(co2) |>
  ggplot(aes(year, co2)) + geom_line() + labs(title = "日本の一人当たりの二酸化炭素排出量（トン）：NA を除く")

df_co2fos |> filter(country == "World") |> drop_na(fossil) |>
  ggplot(aes(year, fossil)) + geom_line() + labs(title = "世界の化石燃料による発電量（％）")

化石燃料のデータは、1970年からデータがあることも分かりましたね。データの数を最初に調べたときのことを思い出してください。特に、色付きの方。違う条件で描くと、比較ができます。いろいろなことが分かりますよ。

ここでは、日本と世界と変えましたが、同じ国または地域で比較するのが適切でしょう。いくつもの国や、地域でみてみることはおすすめです。

5.2.1.1.1 参考

では、次のようにしたらどうでしょう。

df_co2fos |> ggplot(aes(year, fossil)) + geom_line() +
  labs(title = "Filter を使わなかった時の折線グラフ")

なにか変ですね。これは、それぞれの年にたくさんデータがあって、それを線では結べないので、こうなってしまうのです。実は、geom_line を geom_point に変えれば、多少状況は分かります。

df_co2fos |> ggplot(aes(year, fossil)) + geom_point() +
  labs(title = "Filter を使わなかった時の散布図")

これを、折線で結ぶことはできません。また、missing value もたくさんあるよとなっています。そこで、ちょっと改良します。

df_co2fos |> filter(!(iso2c %in% REGION)) |> drop_na(fossil, income) |> 
  ggplot(aes(year, fossil, col = factor(income, levels = INCOME))) + geom_point() +
  labs(title = "国別、所得レベル別の一人当たりの化石燃料での発電割合", col = "Income Level")

df_co2fos |> filter(!(iso2c %in% REGION)) |> drop_na(co2, income) |> 
  ggplot(aes(year, co2, col = factor(income, levels = INCOME))) + geom_point() +
  labs(title = "国別、所得レベル別の一人当たりの二酸化炭素排出量", col = "Income Level")

現時点では、何をやっているかわからないかもしれませんが、二つのグラフは、散布図のひとつで、所得レベルで、色を変えています。また、その表示の順序を、何も指定しないと、アルファベット順になるので、ここでは、特別に、INCOME に登録した順序で表示しています。いろいろなことを見てとることもできると思います。あなたは、どんなことに気づきましたか。

実は、散布図の方が基本的です。そちらを先に学ぶ方がよいように思います。ただ、世界開発指標は、年毎のデータになっていますから、経年変化を先にしています。上に書いたように、経験変化を調べる時は、まずは、国などを指定する必要があることを、理解してください。

5.2.1.2 国または地域の二つ以上の変数に関する経年変化

このためには、縦長の表を使います。

ただ、この例では、二つの指標の単位がことなるため、このように表示するのが適切かどうかは不明です。

表題が長い時は、‘’ を挿入すれば、改行できます。

df_co2fos_long |> filter(country == "Japan") |> drop_na(value) |>
  ggplot(aes(year, value, col = name)) + geom_line() + 
  labs(title = "日本の一人当たりの二酸化炭素排出量と、\n化石燃料による発電の割合")

df_co2fos_long |> filter(country == "World") |> drop_na(value) |>
  ggplot(aes(year, value, col = name)) + geom_line() +
  labs(title = "世界の一人当たりの二酸化炭素排出量と、\n化石燃料による発電の割合")

5.2.1.2.1 参考

実は、指標の数が二つぐらいであれば、df_co2fos を使って、描くことも可能です。座標軸名を変更しています。

df_co2fos |> filter(country == "Japan") |> drop_na(co2, fossil) |>
  ggplot(aes(x = year)) + geom_line(aes(y = co2), col = "red") + 
  geom_line(aes(y = fossil), col = "blue") +
  labs(title = "df_co2fos を使って一つ一つの指標を描画", y = "fossil (%) in blue, co2 (ton) in red")

co2 の方は、大体、9ton ぐらいなので、7倍ぐらいしてみましょう。

df_co2fos |> filter(country == "Japan") |> drop_na(co2, fossil) |>
  ggplot(aes(x = year)) + geom_line(aes(y = co2*7), col = "red") + geom_line(aes(y = fossil), col = "blue") + labs(title = "df_co2fos を使って一つ一つの指標を描画：co2 七倍", y = "fossil (%) in blue, co2 (ton) x 7 in red")

どうですか、縦軸は、少し変ですが、いままでとは、違って気づくこともありますよね。両方のデータがあるところだけを表示させています。

簡単なのは、df_co2fos_long を使った方でしたが、df_co2fos を使った方も、使えそうなことが分かりました。

色の取り扱いもいろいろとありますが、それは、参考文献などを参照してください。

5.2.1.3 データの確認

グラフが表示されなかったり、期待したものではない時はそのデータを見てみましょう。 途中から、co2 と fossil の両方のデータがあることがわかります。上の方で、各年のデータの数の棒グラフでも一般的傾向を確認してありますね。

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

5.2.1.4 いくつかの国または地域の二つ以上の変数に関する経年変化

regions は、国名で定義していましたから、filter(country %in% regions) で選択しています。次の、ASEAN は、iso2c で定義していましたから、filter(iso2c %in% ASEAN) としています。違いを理解してください。

filter の使い方、一つの国や地域のときは、filter(country == “Japan”), filter(country == “World”) でした。今度は、%in% をつかい、そのあとのグループに入っているものを選択しています。

参考：

regions <- c("South Asia", "Europe & Central Asia", "Middle East & North Africa", 
"East Asia & Pacific", "Sub-Saharan Africa", "Latin America & Caribbean", 
"North America")

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

df_co2fos |> filter(country %in% regions) |> drop_na(co2) |>
  ggplot(aes(year, co2, col = country)) + geom_line() +
  labs(title = "地域ごとの一人当たりの二酸化炭素排出量")

選択は、iso2c でしていますが、col = country としてあるので、国名が表示されています。この次の 一人当たりの二酸化炭素排出量の方のグラフでは、それを、iso2c に変更しています。

df_co2fos |> filter(iso2c %in% ASEAN) |> drop_na(fossil) |>
  ggplot(aes(year, fossil, col = country)) + geom_line()  +
  labs(title = "ASEAN の化石燃料による発電割合")

グラフのスペースが確保できる利点もあります。

df_co2fos |> filter(iso2c %in% ASEAN) |> drop_na(co2) |>
  ggplot(aes(year, co2, col = iso2c)) + geom_line()  +
  labs(title = "ASEAN（iso2c 表示）の一人当たりの二酸化炭素排出量")

5.2.1.5 いくつかの国または地域の二つ以上の変数に関する経年変化

上にも書いたように、単位が異なるので、適切かは疑問ですが、同時にみることはできます。

df_co2fos_long |> filter(country %in% regions) |> drop_na(value) |>
  ggplot(aes(year, value, col = country, linetype = name)) + geom_line() +
  labs(title = "一人当たりの二酸化炭素排出量と\n化石燃料による発電割合")

df_co2fos_long |> filter(country %in% SACU) |> drop_na(value) |>
  ggplot(aes(year, value, col = country, linetype = name)) + geom_line() +
  labs(title = "南部アフリカ関税同盟の一人当たりの二酸化炭素排出量と\n化石燃料による発電割合")

5.3 分布

それぞれの指標の分布をみます。比較するには、箱ひげ図も有効です。

5.3.1 ヒストグラム

bins 区分けの数、binwidth 区分幅を変更して、適切に表示されるようにします。何も書いていないと、初期設定では、bins = 30 になっています。

5.3.1.1 例

特定の年に、一定の範囲のところにはいる、国の数（度数分布）を表します。

df_co2fos |> filter(!(iso2c %in% REGION)) |> filter(year == 2020) |> drop_na(co2) |>
  ggplot(aes(co2)) + geom_histogram() +
  labs(title = "2020年の国別の一人当たりの二酸化炭素排出量のヒストグラム")

区分けの個数を10個 bins = 10 に変更すると次のようになります。

df_co2fos |> filter(!(iso2c %in% REGION)) |> filter(year == 2020) |> drop_na(co2) |>
  ggplot(aes(co2)) + geom_histogram(bins = 10) +
  labs(title = "2020年の国別の一人当たりの二酸化炭素排出量のヒストグラム")

区分けの幅を binwidth = 2 とすると下のようなグラフになります。

df_co2fos |> filter(!(iso2c %in% REGION)) |> filter(year == 2020) |> drop_na(co2) |>
  ggplot(aes(co2)) + geom_histogram(binwidth = 2) +
  labs(title = "2020年の国別の一人当たりの二酸化炭素排出量のヒストグラム")

地域名を色で加えることも可能です。塗りつぶしなので、fill を使っています。

df_co2fos |> filter(!(iso2c %in% REGION)) |> filter(year == 2020) |> drop_na(co2) |>
  ggplot(aes(co2, fill = region)) + geom_histogram(bins = 10) +
  labs(title = "2020年の国別の一人当たりの二酸化炭素排出量のヒストグラム\n地域情報を加えたもの")

境目がみにくいように思うので、わたしは、黒枠の細めの線を入れます。

df_co2fos |> filter(!(iso2c %in% REGION)) |> filter(year == 2020) |> drop_na(co2) |>
  ggplot(aes(co2, fill = region)) + geom_histogram(bins = 10, col = "black", linewidth = 0.2)+
  labs(title = "2020年の国別の一人当たりの二酸化炭素排出量のヒストグラム\n地域情報を加えたもの")

左に偏っているので、LOG 表示を使うこともひとつです。

df_co2fos |> filter(!(iso2c %in% REGION)) |> filter(year == 2020) |> drop_na(co2) |>
  ggplot(aes(co2, fill = region)) + 
  geom_histogram(bins = 10, col = "black", linewidth = 0.2) + scale_x_log10()+
  labs(title = "2020年の国別の一人当たりの二酸化炭素排出量のヒストグラム\n地域情報を加えたもの", x = "co2 (ton/person, log10 scale)")

もう一つの変数、fossil についてもやってみてください。

5.3.1.1.1 参考

年を、2020 にしましたが、fossil のときにも同じように、2020 にするとどうなるでしょうか。

df_co2fos |> filter(!(iso2c %in% REGION)) |> filter(year == 2020) |> drop_na(fossil) |>
  ggplot(aes(fossil)) + geom_histogram(bins = 10) +
  labs(title = "2020年の化石燃料による発電割合")

何も現れません。何が問題なのでしょう。一番最初に、データの数をみましたよね。特に、カラーの方を思い出してください。確かめてみましょう。

df_co2fos |> filter(!(iso2c %in% REGION)) |> filter(year == 2020) |> drop_na(fossil)

2020にすると、NA でないデータはないことが確認できます。では、何年なら大丈夫でしょうか。

df_co2fos |> filter(!(iso2c %in% REGION)) |> filter(year == 2015) |> drop_na(fossil) |>
  ggplot(aes(fossil)) + geom_histogram(bins = 10) +
  labs(title = "2015年の化石燃料による発電割合")

グラフは現れましたが、bins や、binwidth はこれが適切でしょうか。指標によって、変えないといけませんね。

df_co2fos |> filter(!(iso2c %in% REGION)) |> filter(year == 2015) |> drop_na(fossil) |>
  ggplot(aes(fossil, fill = income)) + geom_histogram(binwidth = 10, col = "black", linewidth = 0.2) +
  labs(title = "2015年の化石燃料による発電割合, 所得レベルを含む")

どんなことがわかりますか。

5.3.2 棒グラフ

カテゴリー変数（日本語では名目変数ということばも使われます）と、その値（numerical value）が与えられているときに使います。データの数の時に使った、geom_bar も似ていますが、こちらは、カテゴリー変数や飛び飛びの値（discrete numerical variable）の数値変数などの、数を数えます。

5.3.2.1 例

df_co2fos |> filter(year == 2015, country %in% regions) |> drop_na(co2) |>
  ggplot(aes(country, co2)) + geom_col() +
  labs(title ="地域別一人当たりの二酸化炭素排出量", x = "")

重なり合ってしまうので、軸を変えます。

df_co2fos |> filter(year == 2015, country %in% regions) |> drop_na(co2) |>
  ggplot(aes(country, co2)) + geom_col() + coord_flip() +
  labs(title ="地域別一人当たりの二酸化炭素排出量", x = "")

化石燃料での発電割合はどうでしょうか。次のグラフでは、情報は加えていませんが、塗る色を変えてあります。でも、凡例は不要ですよね。theme(legend.position = “none”) とします。

country は、regions とする手もありますが、消しておいても良さそうです。注意が必要なのは、coord_flip をする前の座標なので、x を消します。theme(legend.position = ’none”)

df_co2fos |> filter(year == 2015, country %in% regions) |> drop_na(fossil) |>
  ggplot(aes(country, fossil, fill = country)) + geom_col() + coord_flip()

df_co2fos |> filter(year == 2015, country %in% regions) |> drop_na(fossil) |>
  ggplot(aes(country, fossil, fill = country)) + geom_col() + coord_flip()  +
  labs(title ="地域別化石燃料による発電量（%）", x = "", y = "") +
  theme(legend.position = "none")

fill（塗りつぶし）を使いましたが、col だとどうなるのでしょうか。あまり適切だとは思えません。しかし、ヒストグラムのときのように、枠を黒にすることは、多少意味があるかもしれません。

df_co2fos |> filter(year == 2015, country %in% regions) |> drop_na(fossil) |>
  ggplot(aes(country, fossil, col = country)) + geom_col() + coord_flip() + 
  labs(col = "")

df_co2fos |> filter(year == 2015, country %in% regions) |> drop_na(fossil) |>
  ggplot(aes(country, fossil, fill = country)) + geom_col(col = "black", linewidth = 0.2) +
  coord_flip() + labs(title = "地域別化石燃料による発電量（%）", x = "", y = "") + theme(legend.position = "none")

5.3.2.1.1 参考

実は、縦のままにすることも可能です。ただ、それには、strngr という、tidyverse の一部ではありますが、読み込まれていない、パッケージの命令を使います。

df_co2fos |> filter(year == 2015, country %in% regions) |> drop_na(fossil) |>
  ggplot(aes(country, fossil, fill = country)) + geom_col(col = "black", linewidth = 0.2) +
  scale_x_discrete(labels = function(x) stringr::str_wrap(x, width = 15)) +
  labs(title = "地域別化石燃料による発電量（%）", x = "", y = "") + theme(legend.position = "none")

5.3.2.2 データが十分ある最近の年の値の10カ国の値の棒グラフ

表示の順序を制御するのには、factor というものを使います。カテゴリー変数に順序をつけるものです。普通は、factor(var, legends = c(var の 順序)) とします。ただ、基本的な順序付を制御することは、個別に定義しなくても可能です。それが、fct_inorder で、データに現れる順序に表示します。ここでは、大きい順に表示されます。それを、fct_rev で逆にしています。fct_rev をとった、この次の例を比較してください。

df_co2fos |> filter(!(iso2c %in% REGION)) |> filter(year == 2015) |> drop_na(co2) |>
  arrange(desc(co2)) |> slice_head(n = 10) |>
  ggplot(aes(fct_rev(fct_inorder(country)), co2)) + geom_col() + coord_flip() + labs(title = "一人当たりの二酸化炭素排出量が多い順", x = "")

df_co2fos |> filter(!(iso2c %in% REGION)) |> filter(year == 2015) |> drop_na(co2) |>
  arrange(desc(co2)) |> slice_head(n = 10) |>
  ggplot(aes(fct_inorder(country), co2)) + geom_col() + coord_flip() + labs(title = "一人当たりの二酸化炭素排出量が多い順", x = "")

地域ことに色分けすることも可能です。

df_co2fos |> filter(!(iso2c %in% REGION)) |> filter(year == 2015) |> drop_na(co2) |>
  arrange(desc(co2)) |> slice_head(n = 10) |>
  ggplot(aes(fct_rev(fct_inorder(country)), co2, fill = region)) + geom_col() + coord_flip() + labs(title = "一人当たりの二酸化炭素排出量が多い順", x = "")

df_co2fos |> filter(!(iso2c %in% REGION)) |> filter(year == 2015) |> drop_na(fossil) |>
  arrange(fossil) |> slice_head(n = 10) |>
  ggplot(aes(fct_rev(fct_inorder(country)), fossil)) + geom_col() + coord_flip()  + labs(title = "化石燃料による発電割合の少ない順", x = "")

所得レベルで色分けをしてみましょう。

df_co2fos |> filter(!(iso2c %in% REGION)) |> filter(year == 2015) |> drop_na(fossil) |>
  arrange(fossil) |> slice_head(n = 10) |>
  ggplot(aes(fct_rev(fct_inorder(country)), fossil, fill = income)) + geom_col() + coord_flip()  + labs(title = "化石燃料による発電割合の少ない順", x = "")

5.3.3 箱ひげ図（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()

5.3.3.1 例

国や地域はカテゴリー変数ですから、地域別で、一人あたりの二酸化炭素排出量を箱ひげ図に書いてみましょう。ここでは、df_co2fos_long を使っていますが、このような箱ひげ図は、df_co2fos でも同様に描けます。ただ、コードは異なりますから注意してください。

df_co2fos_long のときは、name に、co2 と fossil が含まれていますから、まず、co2 を選択します。グラフを書く時は、value にある値を使います。

df_co2fos_long |> filter(name == "co2") |> 
  filter(region != "Aggregates") |> drop_na(value) |>
  ggplot(aes(region, value)) + geom_boxplot() + coord_flip() +
  labs(title = "地域ごとの一人当たりの二酸化炭素排出量（トン）")

今度は、df_co2fos を使います。こちらでは、filter(region != “Aggregates”) としていますが、REGION を使うてもあります。それは、次に示します。実は、次の方が正確です。それは、ちょっと細かい話ですが、地域名が、入っていない国があるのです。そのこともあって、REGION を最初に定義してあります。

df_co2fos |> filter(region != "Aggregates") |> drop_na(fossil) |>
  ggplot(aes(region, fossil)) + geom_boxplot() + coord_flip() +
  labs(title = "地域ごとの化石燃料による発電量の割合")

REGION は、iso2c で定義していましたから、filter(!(iso2c %in% REGION)) とします。

df_co2fos |> filter(!(iso2c %in% REGION)) |> drop_na(fossil, region) |>
  ggplot(aes(region, fossil)) + geom_boxplot() + coord_flip()+
  labs(title = "地域ごとの化石燃料による発電量の割合")

co2 と、fossil は、一方は、トンで、他方は、パーセントで、尺度が違いますから、同時に、箱ひげ図に表すのは、適切とは言えませんが、一応、df_co2fos_long を使って、コードを書いてみます。

df_co2fos_long |> 
  filter(income != "Aggregates") |> drop_na(value) |>
  ggplot(aes(income, value, fill = name)) + geom_boxplot() + coord_flip() +
  labs(title = "所得レベル別の、化石燃料による発電の割合と\n一人当たりの二酸化炭素排出量")

順序を、INCOME <- c("High income", "Upper middle income", "Lower middle income", "Low income") と指定した順序に帰る時は、factor を使います。

df_co2fos_long |> filter(income != "Aggregates", income != "Not classified") |> 
  drop_na(value) |>
  ggplot(aes(factor(income, levels = INCOME), value, fill = name)) + geom_boxplot() +
  coord_flip() + labs(title = "所得レベル別の、化石燃料による発電の割合と\n一人当たりの二酸化炭素排出量", x = "")

5.4 相関

5.4.1 相関係数

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

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

あくまでも、相関であり、因果関係ではない。再生可能エネルギーを増やしたり、化石燃料による発電を減らしたり、森を増やせば、二酸化炭素排出量が減ると、ここから結論できるわけではない。あくまでも、化石燃料の発電が多い国のほうが、二酸化炭素排出量が多い傾向があるというようなことが見て取れることである。相関関係は因果関係にあらず。Correlation is not causation.

5.4.1.1 例

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

df_co2fos |> filter(!(iso2c %in% REGION)) |> # 地域以外（国のみ選択）
  drop_na(co2, fossil) |> 
  select(co2, fossil) |> cor()

             co2    fossil
co2    1.0000000 0.3195902
fossil 0.3195902 1.0000000

これは、弱い正の相関となります。つまり、「化石燃料による発電が多い国はある程度、二酸化炭素排出量が多い傾向にある」と言えます。

5.4.2 散布図

基本形は、次のようになります。

df_co2fos |> filter で選択 |> drop_na(変数1, 変数2) |>
  ggplot(aes(変数1, 変数2)) + geom_point()

5.4.2.1 例

二変数の相関をみるには、df_co2fos を使います。地域は、filter(!(iso2c %in% REGION)) で、省いてあります。

df_co2fos |> filter(!(iso2c %in% REGION)) |> drop_na(co2, fossil) |>
  ggplot(aes(co2, fossil)) + geom_point() + 
  labs(title = "一人当たりの二酸化炭素排出量と化石燃料での発電割合")

回帰直線を加えます。

df_co2fos |> filter(!(iso2c %in% REGION)) |> drop_na(co2, fossil) |>
  ggplot(aes(co2, fossil)) + geom_point() + geom_smooth(formula = 'y~x', method = "lm", se = FALSE) + ylim(0,100) +
  labs(title = "一人当たりの二酸化炭素排出量と化石燃料での発電割合（回帰直線）")

収入レベルで色分けをすることも可能です。

df_co2fos |> filter(!(iso2c %in% REGION)) |> drop_na(co2, fossil) |> 
  ggplot(aes(co2, fossil, col = factor(income, levels = INCOME))) + geom_point() + 
  labs(title = "一人当たりの二酸化炭素排出量と化石燃料での発電割合", col = "Income Level")

0 に近い方に値がかたまっているので、scale_x_log10() で 対数表時をしてみます。すこし広がります。

df_co2fos |> filter(!(iso2c %in% REGION)) |> drop_na(co2, fossil) |>
  ggplot(aes(co2, fossil)) + geom_point() + scale_x_log10() + 
  labs(title = "一人当たりの二酸化炭素排出量と化石燃料での発電割合")

こちらに回帰直線を加えると下のようになります。

df_co2fos |> filter(!(iso2c %in% REGION)) |> drop_na(co2, fossil) |>
  ggplot(aes(co2, fossil)) + geom_point() + scale_x_log10() + 
  geom_smooth(formula = 'y~x', method = "lm", se = FALSE) + 
  labs(title = "一人当たりの二酸化炭素排出量（LOG）と化石燃料での発電割合")

5.4.3 補足：散布図

一人当たりの購買力平価による国内総生産 GDP (PPP) と、一人当たりの二酸化炭素排出量と、人口のデータを読み込みます。

df_co2_gdp <- WDI(indicator = c(pop = "SP.POP.TOTL", gdppcap = "NY.GDP.PCAP.PP.KD", co2 = "EN.ATM.CO2E.PC"), extra = TRUE)
write_csv(df_co2_gdp, "data/co2_gdp.csv")

df_co2_gdp <- read_csv("data/co2_gdp.csv")

Rows: 16758 Columns: 15── Column specification ────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (7): country, iso2c, iso3c, region, capital, income, lending
dbl  (6): year, pop, gdppcap, co2, 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_co2_gdp

df_co2_gdp |> filter(!(iso2c %in% REGION)) |> drop_na(gdppcap, co2) |> 
  ggplot(aes(gdppcap, co2)) + geom_point() +
  labs(title = "購買力平価による一人当たりの国内総生産と\n一人当たりの二酸化炭素排出量")

明らかに、0 の付近に密集していますから、Log10 スケールを使います。

df_co2_gdp |> filter(!(iso2c %in% REGION)) |> drop_na(gdppcap, co2) |> 
  filter(gdppcap >0, co2 >0) |>
  ggplot(aes(gdppcap, co2)) + geom_point() + scale_x_log10() + scale_y_log10()+
  labs(title = "購買力平価による一人当たりの国内総生産と\n一人当たりの二酸化炭素排出量（ともに Log10 スケール）")

強い相関がみられそうです。回帰直線を引いてみます。subtitle を使ってみました。

df_co2_gdp |> filter(!(iso2c %in% REGION)) |> drop_na(gdppcap, co2) |> 
  filter(gdppcap >0, co2 >0) |>
  ggplot(aes(gdppcap, co2)) + geom_point() + scale_x_log10() + scale_y_log10() +
  geom_smooth(formula = 'y~x', method = "lm", se = FALSE)+
  labs(title = "購買力平価による一人当たりの国内総生産と一人当たりの二酸化炭素排出量", subtitle = "ともに Log10 スケール、回帰直線を加えたもの")

少し点を減らすため、年を、2020 に制限します。さらに、地域を色で加えます。

df_co2_gdp |> filter(!(iso2c %in% REGION), year == 2020) |> drop_na(gdppcap, co2) |> 
  filter(gdppcap >0, co2 >0) |>
  ggplot(aes(gdppcap, co2, col = region)) + geom_point() + scale_x_log10() + scale_y_log10()+
  labs(title = "購買力平価による一人当たりの国内総生産と一人当たりの二酸化炭素排出量",
       subtitle = "ともに 2020年の Log10 スケール のデータ")

さらに、人口の大きさにより、点の大きさを変えます。

df_co2_gdp |> filter(!(iso2c %in% REGION), year == 2020) |> drop_na(gdppcap, co2) |> 
  filter(gdppcap >0, co2 >0) |>
  ggplot(aes(gdppcap, co2, col = region, size = pop)) + geom_point() + scale_x_log10() + scale_y_log10()+
  labs(title = "購買力平価による一人当たりの国内総生産と一人当たりの二酸化炭素排出量",
       subtitle = "ともに 2020年の Log10 スケール のデータ、点の大きさは人口")

5.4.3.0.1 参考

plotly という別のパッケージを使うと、ポインターを点にあてると、情報が読み取れるようになります。

library(plotly)

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

ggp <- df_co2_gdp |> filter(!(iso2c %in% REGION), year == 2020) |> drop_na(gdppcap, co2) |> 
  filter(gdppcap >0, co2 >0) |>
  ggplot(aes(text = country, gdppcap, co2, col = region, size = pop)) + geom_point() + scale_x_log10() + scale_y_log10() +
  labs(title = "一人当たりの国内総生産と二酸化炭素排出量(ポインタでデータ表示)")
ggplotly(ggp)

回帰直線を加えます。

df_co2_gdp |> filter(!(iso2c %in% REGION), year == 2020) |> drop_na(gdppcap, co2) |> 
  filter(gdppcap >0, co2 >0) |>
  ggplot(aes(gdppcap, co2)) + geom_point(aes(col = region, size = pop)) + scale_x_log10() + scale_y_log10() +
  geom_smooth(formula = 'y~x', method = "lm", se = FALSE)+
  labs(title = "購買力平価による一人当たりの国内総生産と一人当たりの二酸化炭素排出量",
       subtitle = "ともに 2020年の Log10 スケール のデータ、点の大きさは人口、回帰直線を含む")

5.4.3.0.2 参考：線形モデル概要

直線とフィットしているかをみるのは、線形回帰分析と言いますが、その概要は以下のようになります。

上は、対数をとる前と、下は、対数をとってから。

第６週に、少し紹介しましたから、大体どのようなことを言っているかだけ、書いておきます。

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

df_co2_gdp |> filter(!(iso2c %in% REGION), year == 2020) |> drop_na(gdppcap, co2) |> 
  filter(gdppcap >0, co2 >0) |> lm(co2 ~ gdppcap, data = _) |> summary()

Call:
lm(formula = co2 ~ gdppcap, data = filter(drop_na(filter(df_co2_gdp, 
    !(iso2c %in% REGION), year == 2020), gdppcap, co2), gdppcap > 
    0, co2 > 0))

Residuals:
    Min      1Q  Median      3Q     Max 
-9.3224 -1.0320 -0.6497  0.3410 16.1035 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 6.194e-01  3.351e-01   1.849   0.0662 .  
gdppcap     1.686e-04  1.203e-05  14.012   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.251 on 180 degrees of freedom
Multiple R-squared:  0.5217,    Adjusted R-squared:  0.519 
F-statistic: 196.3 on 1 and 180 DF,  p-value: < 2.2e-16

df_co2_gdp |> filter(!(iso2c %in% REGION), year == 2020) |> drop_na(gdppcap, co2) |> 
  filter(gdppcap >0, co2 >0) |> lm(log10(co2) ~ log10(gdppcap), data = _) |> summary()

Call:
lm(formula = log10(co2) ~ log10(gdppcap), data = filter(drop_na(filter(df_co2_gdp, 
    !(iso2c %in% REGION), year == 2020), gdppcap, co2), gdppcap > 
    0, co2 > 0))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.66078 -0.16519 -0.01013  0.16526  0.61827 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)    -4.33682    0.15238  -28.46   <2e-16 ***
log10(gdppcap)  1.13770    0.03739   30.43   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.249 on 180 degrees of freedom
Multiple R-squared:  0.8373,    Adjusted R-squared:  0.8364 
F-statistic:   926 on 1 and 180 DF,  p-value: < 2.2e-16

この線形回帰モデルによる分析は、いくつかの仮定の上に載っていますので、個人的には、現時点では利用をお勧めしません。ただ、これから、分野によっては、この情報をよく使うことになると思いますので、大雑把なことだけ、上に紹介しました。

df_co2_gdp |> filter(!(iso2c %in% REGION), year == 2020) |> drop_na(gdppcap, co2) |> 
 select(gdppcap, co2) |> cor()

          gdppcap       co2
gdppcap 1.0000000 0.7222794
co2     0.7222794 1.0000000

0.7222794*0.7222794

[1] 0.5216875

この値が、log10 をとらない、最初のモデルの、Multiple R Squared の値になっています。

df_co2_gdp |> filter(!(iso2c %in% REGION), year == 2020) |> drop_na(gdppcap, co2) |> 
 filter(gdppcap >0, co2>0) |> select(gdppcap, co2) |> sapply(log10) |> cor()

          gdppcap       co2
gdppcap 1.0000000 0.9150164
co2     0.9150164 1.0000000

0.9150164 * 0.9150164

[1] 0.837255

確かに、この値が、log10 をとった方の、Multiple R Squared の値になっています。

6 コミュニケーション

探索的データ分析の最後は、コミュニケーションです。最後には、プレゼンテーション（発表）や、レポート（paper）の作成になります。わかりやすく伝えるには、作業履歴がしっかり書かれている、R Notebook よりも、簡潔にまとめたものが必要とされます。探索的データ分析は、再現性（Reproducibility）が大切ですから、実際の作業は、R Notebook で行うことを紹介していきました。これによって、根拠が求められば、どのように、データを取得し、どのように、図が表示されたかも、簡単に、シェアすることができます。しかし、分析結果、それから、わかることを正確に、かつ、簡潔に伝えることも必要です。そのために、どうしても必要なことを書いておきます。

一般的な道筋として、

1. スライドで、プレゼンをする。
2. Word または、Google Doc でレポートを書く。
3. PDF にして提出する。

このことを想定して書いていきます。

6.1 図の取り出し

6.1.1 Preview から生成

Preview の図を、Ctrl + Click すると、Copy Image の選択ができます。

Copy Image とすると、図がクリップボードにコピーできると思います。開いていた、Word や、Google Doc に、貼り付けが可能です。サイズもそこで、調節できます。しかし、サイズを変更すると、表題なども一緒にサイズが変わると思います。

6.1.2 Code から生成

1. Code Chunk から、実際のコードの部分をコピーします。（{r} などはコピーしません。）
2. 左下の、Console を開き、ペーストします。そして、実行（Enter など）
3. 右下の、Plots に表示されます。
4. Save as Image … , Save as PDF, Copy to Clipboad と出ます。

• Copy to Clipboard を選択し、Width と、Height などを調節し、Update Preview で確認し、Copy Plot とすれば、Word または、Google doc に貼り付けることができます。
• Save as PDF は、Device Size を調節し、Portrait（縦置き）, Landscape（横置き） を選択し、名前と保存場所を決めて、Save です。
• Save as Image は、画像の種類を選べますが、図のときは、おそらく、png でよいと思います。種類は、圧縮率や圧縮方法が異なります。あとは、サイズと名前と保存場所を決めて、Save です。

保存したものは、Word や、Google Doc から、Image の挿入を選択し、画像を選択することになります。図は一つのときは、良いですが、幾つもあるとわからなくなりますから、適切に、名前をつけることも大切です。

6.2 文書の作成

日本語の場合は特にいくつか注意することがありますので、まとめておきます。Windows, Mac, Linux で、方法が違う部分もありますが、極力、どれでも同じにできる方法を紹介します。（分断をさけるため）

6.2.1 R Markdown テンプレート

• R Markdown いろいろな形式での出力 [リンク], [Rmd]

このファイルを使うと、さまざまな形式に出力できるだけではなく、説明も書かれていますので、まずは、これを丁寧に読んでください。基本だけ、少し書いておきます。

6.2.2 Word への出力

上の knit ボタンから、Word を作成できます。しかし、日本語が、図の表題に入っている時は、文字化けが起きる可能性がありますので、必ず、次を入れておいてください。

library(showtext) 
knitr::opts_chunk$set(fig.showtext=TRUE)

また、初期設定では、紙のサイズが、US Letter になっています。印刷するときなど、問題が起きる可能性がありますから、必ず、A4 に修正してください。このような問題をさけるためにも、基本的に、R Markdown（R Notebook はその一形式）で文章を作り、Word に出力する場合には、形式を整えた、参照ファイルを作って、それにそって、出力するようにすることもできます。詳しくは、[R Markdown クックブック] を参照してください。

R Notebook から、Word を作成し、それを編集することで、十分だと思います。

最終的に、PDF にする時は、Word の機能の、Print to PDF を使って、印刷を選び、PDF に保存をします。

直接、PDF を作成することも可能ですが、さまざまな制御が複雑ですから、お勧めしません。基本的なことは、上に引用した、「R Markdown いろいろな形式での出力」に書いてあります。

6.2.3 スライド作成

基本的な部分は、「R Markdown いろいろな形式での出力」に書いてあります。参照してください。

PowerPoint を使われる方が多いかもしれませんが、そのテンプレートでの制御も可能です。ただ、個人的には、PowerPoint のテンプレート作成は手間がかかるので、これも、PowerPoint に出力して、そちらを修正した方がお勧めです。

わたしは、授業では、isoslides を利用しています。軽いのと、HTML 形式なので、表の、paged 機能も使えるため、見せながら、ページ送りもできるからです。

いずれにしても、スライド作成には、ページサイズにあった、ページ送りの場所をきめないといけません。## の Heading 2 があると、ページ替えができますが、それ以外では、— ハイフン三つ、または、Visual Editor の Insert から、Horizontal Rule を選べば、回ページできます。

ただ、Horizontal Rule の後には、空白行が必要です。それがないと、Error になります。

6.2.4 さまざまな形式

上の、Preview や、knit の右にある、ギアマークの一番下に、Output Options があります。図のサイズや、他の項目も変更可能です。

図は、初期設定では、6.5 in x 4 in となっています。また、それぞれの Code Chunk の右上のギアマークから、それぞれの図のサイズも指定できます。その文書の図のサイズの初期設定をして、あとは、個別の図のサイズを変更という具合です。

6.2.5 コード・チャンク

上に個別の図のサイズの変更は、コードチャンクの右上のギアマークからと書きました。同じ、ギアマークから、いくつかの変更が可能です。簡単に書いておきます。

• Chunk Name: 名前をつけておくと、エラーが起こった時にどこだかすぐわかりますが、同じ名前のチャンクが二つあると、エラーになるので、わたしは、つけていません。
• Output Option: (Use document default), Show output only, show code and output, show nothing (run code), show nothing (don’t run code)
  • (Use document default) 文書の設定通りという意味ですが、特に設定してなければ、関係ありません。
  • Show output only: これは、コードは見せないときに使います。文書を作成する時には、コードはいらないことも多いでしょう。
  • show code and output: これが通常の状態です。初期設定で帰る場合もあるので、このオプションも入れてあるのでしょう。
  • show nothing (run code): コードを表示せず、Warning や、message も出さず、実行結果だけを示すということです。Show output only ではなく、こちらを使った方がよい場合が多いかと思います。ただ、それは、最後にWord などの文書を作成する時で、最初からこれにしてあると、問題がわかりません。
  • show nothing (don’t run code): エラーが修正できず、削除したい場合に選択します。実は、これをすると、{r eval=FALSE, include=FALSE} とはいるのですが、コードはみせるが、実行はしない時は、eval = FALSE だけを入れておくように、わたしはしています。

最後の Chunk Options を選択すると、Option の詳細が書かれています。

7 おわりに

7.1 トラブル・シューティング

いろいろな問題が背景にあり、一般的にはかけません。世界の問題・課題と同じです。これが、問題ですということはできませんよね。

しかし、演習を担当して、質問されたり、課題をみていて、わかった問題などを中心に少し書いておきます。

• Error: … does not exist.
  • さまざまな背景が考えられます。まずは、この最初に説明した、R Notebook の保存場所に、data フォルダーがないために、書き込めないことがあります。大切なのは、R Notebook と、data フォールダがあるかないかではなく、どこにあるか、同じ場所にあるかどうかです。多くの場合、data フォルダの中に、R Notebook が作成されてて、data に書き込めないという状況が生じていました。write_csv(df_co2_fossil, “data/co2_fossil.csv”) は、読み込んだ、df_co2_fossil と名前をつけたデータを、data フォルダの中に、co2_fossil.csv という名前の、CSV 形式（カンマ区切られたテキストデータ）として、書き込むという命令です。data フォルダがなければ、書き込めませんし、data フォルダが、R Notebook のある場所に、なければ書き込めません。次に、df_co2_fossil <- read_csv(“data/co2_fossil.csv”) として、もう一度読み込んでいますが、このときも、data フォルダのなかに、co2_fossil.csv という名前の、CSV ファイルがないと読み込めません。
  • 何かがないと言われたら、それが最初に作成された場所にもどることです。たとえば、df_co2_fossil であれば、最初に、WDI() で読み込むところで、適切に定義されているか、または、 が適切になされているかです。読み込む時に、前のコードのままで、read_csv(“data/co2_fossil.csv”) の部分が、read_csv(“data/unemployment.csv”) などとなっていたら、おかしいですよね。
  • REGION がないと言われたら、それを定義してあるコードが、実行されているか確認する必要があります。そのようなこともあるので、エラーが出たら、Run All Chunks Above をしてみるのも一つです。
  • 次に多いのが、綴りの間違いです。スペルをしっかりチェックしましょう。
• Error: Incomplete expression
  • かっこが足りなかったり、最後に、|> や、+ で終わってる場合によく出ます。|> は次に、データを渡す命令ですから、受け取り側がないわけですし、+ は、ggplot で、図を描く時に、順に繋げていくためのものですから、そのあとには、何かが続くはずです。
• Error: unexpected input in:
  • これは、いろいろな理由があり得ますが、全角空白によって生じる場合があります。Edit から、Find を選び、検索窓に、全角で空白をいれ、Replace の方には、半角で空白を入れ、Next とし、見つかったら、Replace を押すと修正できます。私が例に挙げた、coord_flip() のあとに、全角空白があり、それを、コピーしてしまい、全角空白が入ってしまう例もあったようです。日本語システムでは、無視する場合もありますが、PositCloud ではエラーになります。
• R Notebook で、Run Code をすると、問題なく図が現れるのに、最後に、R Notebook を、Browser で確認すると、最後の方の図が表示されない。
  • これは、途中に、または、--- などが入り込んでいる場合が多いです。わたしは、どうするかというと、Run All をします。すると、エラー箇所がわかる場合があります。そこで、コードを修正する。わからないときは、それより下を、Cut（Windows では、Ctrl + X、Mac では、Command + X）をして、 R Notebook が適切に生成されるかをみます。できていれば、Ctrl + V Mac では、Command + V または、Ctrl + Z (Mac では、Command + Z) で元に戻し、すこし、先からあとを、カットします。と順に、みていくと、どこが問題かがわかります。これは、おそらく、コピーをするときに、間違って、 まで含めてコピーしてしまうようなことから起こるのでしょう。実際、グラフがあると、間違いがわかありづらいこともあるようです。x を押して、図を、消すと、間違いがわかりやすいことが多いようです。せっかく出てきた図を消すのに勇気がいるのでしょうが。
  • 以前、Code chunk の中でも同じようなことをしました。上のほうから、みていき、大丈夫だと思うところの次から、Cut し、実行します。うまくいっていたら、次の塊までもどして、実行。一ステップ、一ステップ、確認することですから、有効で、一般的でもあります。
• なんだかわからなくなったら。
  • 問題は、自分自身がパニックに陥ることです。冷静さをうしない、よぶんなことをどんどんします。まずは、そこまでを　Save し、次に、Sava as で違った名前にして、カットなどを使って、上から順にみていく。それでも、うまくいかなかったら、最初からやり直す。ただし、大丈夫そうなところは、コピー・ペーストして、前のものを使う。
  • 編集にもある、Undo, Redo も有効ですよ。だいたい、間違った時は、直近で変なことをしたことが多いので。

（また気づいたら書いておきます）

7.2 エラー対応

• 入力したときには、例を参照して、スペルを確認。全角や変数名の空白は問題を起こす。括弧　() や引用符がペアでマッチしているか、確認。引用符が入っていなかったり、== のところが、= だったり、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 は、かならず残すこと。他にも、select で短いデータにしたり、long データを作成するようなときに、{r} の部分を、{r cache = TRUE} として、キャッシュすることで、メモリーを節約することで、エラーを回避できる場合もあります。しかし、おすすめなのは、前半と、後半などに分けて、ファイルを作成することです。最初のデータをよみ込むところは必要ですから、まず、コピーを作成し、最初のファイルは視覚化の後半を消去、二つ目のファイルは、視覚化の前半を消去とするなどしてみてください。短いファイルであれば、問題ないと思います。

7.3 参考文献

• 課題1：解説 リンク・解答例 リンク
• 第6週演習：解答例 リンク

1. 「みんなのデータサイエンス - Data Science for All」[はじめてのデータサイエンス]

2. Posit Recipes（旧 Posit Primers）: The Basics 対話型の演習サイトの最初 [Link]

3. Posit Cheat Sheet. 早見表です。印刷して使うために、PDF も提供しています。[Site Link]

4. DataCamp Cheat Sheet: Tidyverse for Biginners. データサイエンスの教育をしている会社の早見表の一つです。基本が簡単にまとまっています。[Link]

5. 箱ひげ図の見方 Link、外れ値検出のある箱ひげ図 Link

6. Video: Histogram, Boxplot