ECO232 Rではじめるデータ・サイエンス

1 DAY 1：2月20日

1.1 RStudio で Project の作成

1.1.1 確認

1. R と R Studio のインストール
2. RStudio: New Project 作成

• R Studio の起動
• New Project: eco232, eco232.Rproj

3. Save して、終了
4. プロジェクトを起動確認

1.2 基本コマンド - in Console

• head(cars)
• str(cars)
• summary(cars)
• df_cars <- cars
  • <- と _ と %>% と ` back tick の確認
• View(cars) または、右上の Environment から、df_cars をクリック
• ?cars または Help 検索窓で cars, head など

おすすめ：Sys.setenv(LANG = "en")

1.3 Package の確認とインストール

R packages are extensions to the R statistical programming language containing code, data, and documentation in a standardised collection format that can be installed by users of R using Tool > Install Packages in the top menu bar of R Studio.

Rパッケージは、Rの拡張機能で、コード、データ、ドキュメントを標準化されたコレクション形式で含んでおり、標準的なものは、R Studio の Top Bar の Tool > Install Packages からインストールできます。

• Minimum: tidyverse, rmarkdown, WDI

あとから使うので、ロードしておきます。最初に次のようなコードを実行します。右の三角を押します。

library(tidyverse)
library(WDI)

1.4 R Markdown 入門

1.4.1 R Notebook

R Markdownはデータサイエンスのためのオーサリングフレームワーク。

コード（プログラム）とその実行結果、を記録・表示し、高品質のレポートの作成を可能にします。

R Notebook は、独立してインタラクティブに実行できるチャンクを持つR Markdownドキュメントの一つの形式で、入力のすぐ下に出力が表示することができます。

1. File > New File > R Notebook
2. Save with a file name, say, test-notebook
3. Preview by [Preview] button
4. Run Code Chunk plot(cars) and then Preview again.

1.5 World Development Indicator (WDI)

WDI(country = "all", indicator = c(gdp = "NY.GDP.MKTP.CD"),
    extra=TRUE) %>% drop_na(gdp) %>%
  filter(year==max(year), income !="Aggregates") %>% 
  drop_na(region) %>% arrange(desc(gdp))

chosen_countries <- c("United States","China", "Japan", "Germany", "United Kingdom","India")
WDI(country = c("CN","GB","JP","IN","US","DE"), indicator = c(gdp = "NY.GDP.MKTP.CD"), extra=TRUE) %>% drop_na(gdp) %>% 
  ggplot(aes(year, gdp, col = country)) + geom_line() +
  labs(title = "WDI NY.GDP.MKTP.CD: gdp")

WDI(country = c("CN","IN","JP","US"), 
    indicator = c(gdp_growth_rate = "NY.GDP.MKTP.KD.ZG"), extra=TRUE) %>%
  drop_na(gdp_growth_rate) %>% 
  ggplot(aes(year, gdp_growth_rate, col = country)) + geom_line() +
  labs(title = paste("WDI NY.GDP.MKTP.KD.ZG: gdp growth rate"))

1.5.1 指標 Indicators (WDI)

• NY.GDP.MKTP.CD: GDP (current US$)
• NY.GDP.MKTP.KD.ZG: GDP growth (annual %)

1.5.2 指標 WDI を探してみよう

The World Development Indicators is a compilation of relevant, high-quality, and internationally comparable statistics about global development and the fight against poverty. The database contains 1,400 time series indicators for 217 economies and more than 40 country groups, with data for many indicators going back more than 50 years.

WDIは、世界の開発状況と、貧困との戦いに関する、適切で上質、かつ、国際的に比較可能な時系列の統計データを編纂したものです。このデータベースは、217の経済と40以上の国グループについて1,400の時系列指標を含み、指標のデータの多くは50年以上前に遡ることができます。

• 世界銀行（World Bank）:
• World Bank Open Data:
  • Country / Indicator > Featured & All > Details
• World Development Indicators (WDI) :
  • Themes: Poverty and Inequality, People, Environment, Economy, States and Markets, Global Links
  • Open Data & DataBank: Explore data, Query database

1.5.3 指標 WDIの例

• NY.GDP.MKTP.CD: GDP (current US$)
• NY.GDP.DEFL.KD.ZG: Inflation, GDP deflator (annual %)
• SL.UEM.TOTL.NE.ZS: Unemployment, total (% of total labor force) (national estimate)
• CPTOTNSXN: CPI Price, nominal
• SL.TLF.CACT.MA.NE.ZS: Labor force participation rate, male (% of male population ages 15+) (national estimate)
• SL.TLF.CACT.FE.NE.ZS: Labor force participation rate, female (% of male population ages 15+) (national estimate)

1.5.4 練習 1. - 調べてみたい WDI

いくつか、リストしてみましょう。

1.6 WDI パッケージ

WDI パッケージで、データをダウンロードしたり、探したり、詳細情報を得たりできます。

1.6.1 指標 WDI 検索

1.6.1.1 検索例 1（WDI名）

WDIsearch(string = "gdp", field = "name", short = TRUE, cache = NULL)

1.6.1.2 検索例 2（WDI）

WDIsearch(string = "NY.GDP.MKTP.CD", field = "indicator", short = TRUE, cache = NULL)

1.6.1.3 練習 2. - 検索（short）

名前で検索（“” の間に、（なるべく簡単な）検索文字列を入れてください。）

WDIsearch(string = "", field = "name", short = TRUE, cache = NULL)

Indicator で検索（“” の間に、調べたい indicator を入れてください。）

WDIsearch(string = "", field = "indicator", short = TRUE, cache = NULL)

1.6.1.4 詳しい情報を得るには

short = FALSE とします。時間がかかるので、検索は、Indicator と、名前などの情報をもったファイルを手元に持っておくことにします。

wdi_cache <- WDIcache()

右上の窓枠（pane）から、wdi_cache を探して、中身を見てみましょう。series と、country の二つのデータ・フレームからなっているリストです。三角印や、右から二番目の巻物のようなアイコンをクリックすると中身が見えます。

1.6.1.5 検索例 3（WDI名）

WDIsearch(string = "CPI Price", field = "name", short = FALSE, cache = wdi_cache)

• CPTOTNSXN: CPI Price, nominal
  • The consumer price index reflects the change in prices for the average consumer of a constant basket of consumer goods. Data is not seasonally adjusted.

1.6.1.6 検索例 4（WDI）

WDIsearch(string = "NY.GDP.MKTP.KD.ZG", field = "indicator", short = FALSE, cache = wdi_cache)

1.6.1.7 練習 2 - 検索（long w/ cache）

string と、field を、ふたつとも入れてください。

WDIsearch(string = "", field = "", short = FALSE, cache = wdi_cache)

2 DAY 2：2月22日

2.1 （続）WDI パッケージ

2.1.1 指標 WDI データのダウンロード

Indicator が決まったら、ダウンロードします。

?WDI

2.1.1.1 ダウンロード例 1-1

df_gdp1 <- WDI(country = "all", indicator = "NY.GDP.MKTP.CD")
df_gdp1

2.1.1.2 ダウンロード例 1-2

df_gdp2 <- WDI(country = "all", indicator = c(gdp = "NY.GDP.MKTP.CD"))
df_gdp2

2.1.1.3 ダウンロード例 1-3

df_gdp3 <- WDI(country = "all", indicator = c(gdp = "NY.GDP.MKTP.CD"), extra=TRUE, cache=wdi_cache)
df_gdp3

2.1.1.4 ダウンロード例 1-4

df_gdp4 <- WDI(country = c("CN","GB","JP","IN","US","DE"), indicator = c(gdp = "NY.GDP.MKTP.CD"), extra=TRUE, cache=wdi_cache)
df_gdp4

2.1.1.5 ダウンロード例 2-1

• NY.GDP.DEFL.KD.ZG: Inflation, GDP deflator (annual %)
• CPTOTNSXN: CPI Price, nominal

df_gdp21 <- WDI(country = "all", 
                indicator = c(gdp_deflator = "NY.GDP.DEFL.KD.ZG", 
                              cpi_price = "CPTOTNSXN"), 
                extra=TRUE, cache=wdi_cache)
df_gdp21

str(df_gdp21)

'data.frame':   23972 obs. of  14 variables:
 $ country     : chr  "Advanced Economies" "Advanced Economies" "Advanced Economies" "Advanced Economies" ...
 $ iso2c       : chr  "AME" "AME" "AME" "AME" ...
 $ iso3c       : chr  "" "" "" "" ...
 $ year        : int  1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 ...
 $ status      : chr  "" "" "" "" ...
 $ lastupdated : chr  "2020-07-27" "2020-07-27" "2020-07-27" "2020-07-27" ...
 $ gdp_deflator: num  NA NA NA NA NA NA NA NA NA NA ...
  ..- attr(*, "label")= chr "Inflation, GDP deflator (annual %)"
 $ cpi_price   : num  58.7 60.5 63 66 69.1 ...
  ..- attr(*, "label")= chr "CPI Price,not seas.adj,,,"
 $ region      : chr  NA NA NA NA ...
 $ capital     : chr  NA NA NA NA ...
 $ longitude   : chr  NA NA NA NA ...
 $ latitude    : chr  NA NA NA NA ...
 $ income      : chr  NA NA NA NA ...
 $ lending     : chr  NA NA NA NA ...

summary(df_gdp21)

   country             iso2c              iso3c          
 Length:23972       Length:23972       Length:23972      
 Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character  
                                                         
                                                         
                                                         
                                                         
      year         status          lastupdated       
 Min.   :1960   Length:23972       Length:23972      
 1st Qu.:1982   Class :character   Class :character  
 Median :1996   Mode  :character   Mode  :character  
 Mean   :1995                                        
 3rd Qu.:2009                                        
 Max.   :2021                                        
                                                     
  gdp_deflator         cpi_price         region         
 Min.   :  -98.704   Min.   :  0.00   Length:23972      
 1st Qu.:    2.317   1st Qu.: 55.95   Class :character  
 Median :    5.273   Median : 83.28   Mode  :character  
 Mean   :   25.308   Mean   : 84.18                     
 3rd Qu.:   10.411   3rd Qu.:108.75                     
 Max.   :26765.858   Max.   :551.25                     
 NA's   :11616       NA's   :18410                      
   capital           longitude           latitude        
 Length:23972       Length:23972       Length:23972      
 Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character  
                                                         
                                                         
                                                         
                                                         
    income            lending         
 Length:23972       Length:23972      
 Class :character   Class :character  
 Mode  :character   Mode  :character  
                                      
                                      
                                      
                                      

右上の窓枠の、Environment も見てみましょう。

2.2 可視化 Visualization

グラフ（Chart）を描いて視覚化しよう

2.2.1 グラフ 1

df_gdp4 %>% ggplot(aes(year, gdp, col=country)) + geom_line()

2.2.2 グラフ 2

df_gdp4 %>% drop_na(gdp) %>% 
  ggplot(aes(year, gdp, col=country)) + geom_line() +
  labs(title = paste("WDI - NY.GDP.MKTP.CD: ", "gdp"))

2.2.3 テンプレート Templates

2.2.3.1 一つの国についての、一つの指標（WDI）と、その略称から、折線グラフを作成

Line Plot with one indicator with abbreviation and one country

chosen_indicator <- "SL.UEM.TOTL.NE.ZS"
short_name <- "unemployment"
chosen_country <- "United States"
WDI(country = "all", indicator = c(short_name = chosen_indicator), extra=TRUE, cache=wdi_cache) %>%
  filter(country == chosen_country) %>% 
  ggplot(aes(year, short_name)) + geom_line() +
  labs(title = paste("WDI ", chosen_indicator, ": ", short_name, " - ", chosen_country),
       y = short_name)

2.2.3.2 一つの国についての、一つの指標（WDI）から、折線グラフを作成

Line Plot with one indicator and one country

chosen_indicator <- "SL.UEM.TOTL.NE.ZS"
chosen_country <- "United States"
WDI(country = "all", indicator = c(chosen_indicator = chosen_indicator), 
    extra=TRUE, cache=wdi_cache) %>%
  filter(country == chosen_country) %>% 
  ggplot(aes(year, chosen_indicator)) + geom_line() +
  labs(title = paste("WDI ", chosen_indicator, " - ", chosen_country), 
       y = chosen_indicator)

2.2.3.3 いくつかの国についての、一つの指標（WDI）と、その略称から、折線グラフを作成

Line Plot with one indicator with abbreviation and several countries

chosen_indicator <- "SL.UEM.TOTL.NE.ZS"
short_name <- "unemployment"
chosen_countries <- c("United States","United Kingdom", "Japan")
WDI(country = "all", indicator = c(short_name = chosen_indicator), extra=TRUE, cache=wdi_cache) %>% drop_na(short_name) %>% 
  filter(country %in% chosen_countries) %>% 
  ggplot(aes(year, short_name, col = country)) + geom_line() +
  labs(title = paste("WDI ", chosen_indicator, ": ", short_name), y = short_name)

2.2.3.4 一つの国についての、二つの指標（WDI）と、その略称から、折線グラフを作成

Line Plot with two indicators with abbreviation and one country

chosen_indicator_1 <- "NY.GDP.DEFL.KD.ZG"
short_name_1 <- "gdp_deflator"
chosen_indicator_2 <- "CPTOTSAXNZGY"
short_name_2 <- "cpi_price"
chosen_country <- "United States"
WDI(country = "all", indicator = c(short_name_1 = chosen_indicator_1, short_name_2 = chosen_indicator_2), extra=TRUE, cache=wdi_cache) %>% 
  filter(country == chosen_country) %>% 
  pivot_longer(c(short_name_1, short_name_2), names_to = "class", values_to = "value") %>% drop_na(value) %>%
  ggplot(aes(year, value, col = class)) + geom_line() +
  labs(title = paste("WDI ", chosen_indicator_1, ": ", short_name_1, "\n", chosen_indicator_2, ": ", short_name_2, " - ", chosen_country)) +
  scale_color_manual(labels = c(short_name_1, short_name_2), values = scales::hue_pal()(2))

chosen_indicator_1 <- "SL.TLF.CACT.MA.NE.ZS"
short_name_1 <- "male"
chosen_indicator_2 <- "SL.TLF.CACT.FE.NE.ZS"
short_name_2 <- "female"
chosen_country <- "United States"
WDI(country = "all", indicator = c(short_name_1 = chosen_indicator_1, short_name_2 = chosen_indicator_2), extra=TRUE, cache=wdi_cache) %>% 
  filter(country == chosen_country) %>% 
  pivot_longer(c(short_name_1, short_name_2), names_to = "class", values_to = "value") %>% drop_na(value) %>%
  ggplot(aes(year, value, col = class)) + geom_line() +
  labs(title = paste("WDI ", chosen_indicator_1, ": ", short_name_1, "\n", chosen_indicator_2, ": ", short_name_2, " - ", chosen_country)) +
  scale_color_manual(labels = c(short_name_1, short_name_2), values = scales::hue_pal()(2))

2.2.3.5 いくつかの国についての、二つの指標（WDI）と、その略称から、折線グラフを作成

Line Plot with two indicators with abbreviation and several countries

chosen_indicator_1 <- "NY.GDP.DEFL.KD.ZG"
short_name_1 <- "gdp_deflator"
chosen_indicator_2 <- "CPTOTSAXNZGY"
short_name_2 <- "cpi_price"
chosen_countries <- c("United States", "France", "Japan")
WDI(country = "all", indicator = c(short_name_1 = chosen_indicator_1, short_name_2 = chosen_indicator_2), extra=TRUE, cache=wdi_cache) %>% 
  filter(country %in% chosen_countries) %>% 
  pivot_longer(c(short_name_1, short_name_2), names_to = "class", values_to = "value") %>% drop_na(value) %>%
  ggplot(aes(year, value, linetype = class, col = country)) + geom_line() +
  labs(title = paste("WDI ", chosen_indicator_1, ": ", short_name_1, "\n", chosen_indicator_2, ": ", short_name_2)) +
  scale_linetype_manual(labels = c(short_name_1, short_name_2), values = c("solid", "dashed"))

chosen_indicator_1 <- "SL.TLF.CACT.MA.NE.ZS"
short_name_1 <- "male"
chosen_indicator_2 <- "SL.TLF.CACT.FE.NE.ZS"
short_name_2 <- "female"
chosen_countries <- c("United States", "France", "Japan")
WDI(country = "all", indicator = c(short_name_1 = chosen_indicator_1, short_name_2 = chosen_indicator_2), extra=TRUE, cache=wdi_cache) %>% 
  filter(country %in% chosen_countries) %>% 
  pivot_longer(c(short_name_1, short_name_2), names_to = "class", values_to = "value") %>% drop_na(value) %>%
  ggplot(aes(year, value, linetype = class, col = country)) + geom_line() +
  labs(title = paste("WDI ", chosen_indicator_1, ": ", short_name_1, "\n", chosen_indicator_2, ": ", short_name_2)) +
  scale_linetype_manual(labels = c(short_name_1, short_name_2), values = c("solid", "dashed"))

2.3 課題　Assignment

上のテンプレートをコピーして、下に貼り付け、指標 indicator と、略称 short_name と、いくつかの国名 chosen_countries を、入れ替えて、試してみてください。

2.4 探索的データ解析　Exploratory Data Analysis (EDA)

2.4.1 探索的データ解析とは？ (Posit Primers)

1. EDAは、データが何を語っているかを理解するための反復的なサイクルです。

2. まず、データに関する問いを作成します。

3. データの可視化、変換、モデリングを行い、答えを探します。

学習したことを活用して、問いを修正したり、新しい問いを考えたりします。そして、このサイクルを繰り返していきます。

EDAはデータ分析において重要な役割を果たします。また、データの品質を保証するために、データの質を確認するために使用することもできます。

R4DS からのイメージ

2.4.2 データの取得・読み込み - Importing Data

スタートは、本来は、データの作成・探索ですが、すでに、分析したいデータはすでにあるとして話を進めます。まずは、data フォルダ（directory）を作成しておくと良い。右下の窓枠の Files タブから、New Folder で作成してもよい。

dir.create("./data")

データの取得・読み込みを、四つの方法に分けて説明します。

1. パッケージの利用

• 例：WDI など。何度も、ダウンロードしなくて良いよいに、書き出しておき、2 を使うとよい。write(df_name, "./data/name.csv")

2. コンピュータ上にある CSV などのテキストファイルを読み込む

• 例：df_name <- read_csv("./data/file_name.csv")

3. インターネット上のデータのアドレス（URL）を使って、CSV などのテキストファイルを読み込む。

• 例：df_name <- read_csv(url_of_a_csv)

4. コンピュータ上にある、Excel ファイルなどのデジタルファイルを読み込む。まず、library(readxl)。

• 例：df_name <- read_excel("./data/file_name.xlsx")

5. サイトからダウンロードして、Project のデータフォルダに移す。または、データのアドレス（URL）がわかっていれば、直接ダウンロード。

• 例：`download.file(url_of_a_data, destfile = “./data/data_name”)

6. クリップボードにコピーして読み込む。

• 例：df_name <- read_delim(clipboard())

2.4.3 WDIcache() の扱い

二つの、ファイルが一つになった、リストであるため、違って命令を使います。

wdi_cache <- WDIcache()
write_rds(wdi_cache, "./data/wdi_cache.RData")

wdi_cache <- read_rds("./data/wdi_cache.RData")

2.4.4 国際機関のデータ International Institutions’ Data

• World Bank: https://data.worldbank.org
• UN Data: https://data.un.org
• OECD: https://data.oecd.org/

url_un_pop <- "https://data.un.org/_Docs/SYB/CSV/SYB65_1_202209_Population,%20Surface%20Area%20and%20Density.csv"
df_un_pop0 <- read_csv(url_un_pop)

New names:Rows: 7874 Columns: 7── Column specification ──────────────────────────────────────────────
Delimiter: ","
chr (7): T02, Population, density and surface area, ...3, ...4, .....
ℹ 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_un_pop0

url_un_pop <- "https://data.un.org/_Docs/SYB/CSV/SYB65_1_202209_Population,%20Surface%20Area%20and%20Density.csv"
df_un_pop <- read_csv(url_un_pop, skip=1)

New names:Rows: 7873 Columns: 7── Column specification ──────────────────────────────────────────────
Delimiter: ","
chr (4): ...2, Series, Footnotes, Source
dbl (2): Region/Country/Area, Year
num (1): Value
ℹ 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_un_pop

df_un_pop %>% distinct(`Region/Country/Area`, `...2`)

df_un_pop %>% filter(`Region/Country/Area` %in% c(2,19,142,150,9), Series == "Population mid-year estimates (millions)") %>%
  ggplot(aes(Year, Value, fill = `...2`)) + geom_area(col="black") +
  labs(title = "Population mid-year estimates (millions) of the World")

3 DAY 3：2月24日

3.1 OECD data

• https://data.oecd.org/

3.1.1 例：労働時間当たりGDP

• 日本の時間当たり生産性はOECD38カ国中27位（日本生産性本部「労働生産性の国際比較」）
  • 労働生産性の国際比較2022
• Productivity statistics
  • Read More: Improving Productivity Measurement Practices
    • Level of GDP per capita and productivity
    • GDP per hour worked

Definition of GDP per hour worked

GDP per hour worked is a measure of labour productivity. It measures how efficiently labour input is combined with other factors of production and used in the production process. Labour input is defined as total hours worked of all persons engaged in production. Labour productivity only partially reflects the productivity of labour in terms of the personal capacities of workers or the intensity of their effort. The ratio between the output measure and the labour input depends to a large degree on the presence and/or use of other inputs (e.g. capital, intermediate inputs, technical, organisational and efficiency change, economies of scale). This indicator is measured in USD (constant prices 2010 and PPPs) and indices.

労働時間当たりGDPは、労働生産性の指標である。これは、労働投入量が他の生産要素と組み合わされ、生産プロセスでどれだけ効率的に利用されたかを測定するものである。労働投入量は、生産に従事するすべての人の総労働時間として定義される。労働生産性は、労働者の個人的能力や努力の強さといった労働の生産性を部分的にしか反映していない。アウトプット指標と労働投入量の比率は、他の投入物（資本、中間投入物、技術・組織・効率の変化、規模の経済など）の存在や利用に大きく左右される。この指標は、米ドル（2010年の恒常価格およびPPP）および指標で測定されています。

最初のステップ

1. サイトから、CSV データをダウンロード。
2. プロジェクト内のデータ格納フォルダに入れます。“./data/”
3. 読み込みます。

上のサイトからデータ（Full indicator data）をダウンロードして、プロジェクトの data フォルダーに入れました。

df_oecd_productivity <- read_csv("./data/DP_LIVE_21022023111712065.csv")

Rows: 3894 Columns: 8── Column specification ──────────────────────────────────────────────
Delimiter: ","
chr (6): LOCATION, INDICATOR, SUBJECT, MEASURE, FREQUENCY, Flag Codes
dbl (2): TIME, Value
ℹ 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_oecd_productivity

データの確認

列名を確認し、列ごとに、データを確認

colnames(df_oecd_productivity)

[1] "LOCATION"   "INDICATOR"  "SUBJECT"    "MEASURE"    "FREQUENCY" 
[6] "TIME"       "Value"      "Flag Codes"

unique() のかっこの中に、df_oecd_productivity$LOCATION を入れたものと同じものを出力します。次々とつづけるときに便利なので、このパイプ %>% をわたしはよく使います。unique(df_oecd_productivity$LOCATION) も試してみてください。

df_oecd_productivity$LOCATION %>% unique()

 [1] "AUS"       "AUT"       "BEL"       "CAN"       "CZE"      
 [6] "DNK"       "FIN"       "FRA"       "DEU"       "GRC"      
[11] "HUN"       "ISL"       "IRL"       "ITA"       "JPN"      
[16] "KOR"       "LUX"       "MEX"       "NLD"       "NZL"      
[21] "NOR"       "POL"       "PRT"       "SVK"       "ESP"      
[26] "SWE"       "CHE"       "TUR"       "GBR"       "USA"      
[31] "CHL"       "EST"       "ISR"       "RUS"       "SVN"      
[36] "OECD"      "EU28"      "G-7"       "LVA"       "LTU"      
[41] "EA19"      "ZAF"       "CRI"       "BGR"       "HRV"      
[46] "ROU"       "EU27_2020" "COL"      

df_oecd_productivity$INDICATOR %>% unique()

[1] "GDPHRWKD"

df_oecd_productivity$SUBJECT %>% unique()

[1] "TOT"

df_oecd_productivity$MEASURE %>% unique()

[1] "USD"     "IDX2015"

df_oecd_productivity$FREQUENCY %>% unique()

[1] "A"

df_oecd_productivity$TIME %>% unique()

 [1] 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982
[14] 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995
[27] 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008
[40] 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021

df_oecd_productivity %>% 
  filter(MEASURE == "USD", TIME == 2021) %>% # filtering rows with conditions
  select(LOCATION, Value) %>% # selecting columns
  arrange(desc(Value)) # ordering the rows in the descenting order of the column Value

いくつかの地域を選んで、折線グラフを書いてみる。

３行目がないとどうなりますか。

df_oecd_productivity %>% 
  filter(LOCATION %in% c("JPN", "OECD", "G-7", "EU28")) %>%
  filter(MEASURE == "USD") %>% # same as above up to this line
  ggplot(aes(x = TIME, y = Value, col = LOCATION)) + geom_line() + 
  labs(title="GDP per hour worked", # adding the title and the subtitle
       subtitle="Total, 2015=100, 2021 or latest available")

3.1.2 練習：教育 Education

成人の教育レベルに関する、次のサイトからデータをとって、調べてみてください。

Adult education level: https://data.oecd.org/eduatt/adult-education-level.htm

ダウンロードしたファイルのファイル名が DP_LIVE_21022023120132654.csv と同じであることを確認してください。

df_oecd_education_level <- read_csv("./data/DP_LIVE_21022023120132654.csv")
df_oecd_education_level

上の労働時間当たりGDPの例に習って、データを調べてみてください。

colnames(df_oecd_education_level)

[1] "LOCATION"   "INDICATOR"  "SUBJECT"    "MEASURE"    "FREQUENCY" 
[6] "TIME"       "Value"      "Flag Codes"

df_oecd_education_level$SUBJECT %>% unique()

[1] "BUPPSRY"      "TRY"          "UPPSRY"       "TRY_MEN"     
[5] "TRY_WOMEN"    "UPPSRY_MEN"   "UPPSRY_WOMEN"

これらは、何を意味しているのでしょうか。サイトの Perspectives などをみてください。

他の列（変数）は、どうでしょうか。コードチャンクを挿入して、調べてみてください。

労働時間当たりGDPの例の真似をして、グラフを書いてみました。

df_oecd_education_level %>% 
  filter(LOCATION %in% c("JPN", "OECD", "G-7", "EU28")) %>%
  ggplot(aes(TIME, Value, linetype = SUBJECT, col = LOCATION)) + geom_line() + 
  labs(title="Adult education level")

なにか現れましたが、国も、SUBJECT も期待したものは現れていません。なにが問題なのでしょう。

気づいてたことを書き出してみましょう。

3.2 Excel Data の読み込み

readxl パッケージを使います。これは、tidyverse をインストールするときに、同時に、インストールされますから、インストールは必要ありませんが、tidyverse グループの主たるパッケージはないので、ロードはされていないので、使うときに、ロードする必要があります。

library(readxl) # need to load though readxl is a part of the tidyverse package and installed

実は、範囲を選択し、コピーをし、クリップボードから読む方法もありますが、二つの大きな理由から推奨しません。丁寧に、方法を記述しないと、再現性に問題があること、列のデータ・タイプなどが、適切に、読み込めない場合が多く、読み込んでからの作業が複雑になる。とっても、小さな、自分で作成したデータの場合には、有効かもしれません。

df_excel_clipboard <- read_delim(clipboard())

3.2.1 例: 世界の出生率（Fertility） UN Data

ホットなトピックですね。日本の問題、そして、自分ごととしてだけでなく、世界の状況を見る視点も持ちましょう。

UN data から、検索窓の上にある、Datamarts を選択し、下のほうにある、World Fertility Data United Nations Population Division (UNPD) の [+] 記号を開きます。右の [i] からは情報が得られます。三つのデータがあります。

• Age-specific fertility rates, Total fertility and Mean age at childbearing
• Annual number of live births and Crude birth rate
• Children ever born

まずは、“Age-specific fertility rates, Total fertility and Mean age at childbearing” データを使うことにします。

1. サイトから、Excel データをダウンロード。
2. プロジェクト内のデータ格納フォルダに入れます。“./data/”
3. 試しに、読み込みます。
4. Files から、Import Dataset を選択し、設定をして、もう一度、読み込む。
5. 列名を修正
6. 日本のデータを確認
7. 年齢ごとの数を比較しやすいように、データを変形
8. グラフを書いてみる。
9. 出生平均年齢の表を降年齢順に見てみる。
10. 密度分布グラフを書いてみる。
11. 出生率の表を昇順に見てみる。
12. 密度分布グラフを書いてみる。

Steps 1-3.

1. サイトから、Excel データをダウンロード。
2. プロジェクト内のデータ格納フォルダに入れます。“./data/” {-}
3. 試しに、読み込みます。サイトから、Excel データをダウンロード。{-}

どんなシートがあるか見てみましょう。読み込みに問題があるときは、データ名からスペースを消し、下のコードも修正して読んでみてください。スペースなどがあると問題が起こる場合があります。

excel_sheets("./data/Age-specific fertility rates, Total fertility and .xls")

[1] "FERTILITY INDICATORS" "GENERAL NOTES"       
[3] "SOURCES"             

シートを指定するときは、sheet = 1 とか、sheet = "CHILDREN_EVER_BORN" とします。 指定がないときは、最初のシート。詳細は、Help で、read_excel として、確認してください。たくさんの、オプションがあります。

df_un_fertility <- read_excel("./data/Age-specific fertility rates, Total fertility and .xls")

New names:

df_un_fertility

構造が複雑そうです。

Step 4. Files から、Import Dataset を選択し、設定をして、もう一度、読み込む。

ファイル名も変更しておいが方がよい場合もあります。そのときは、変更記録を RNotebookに追記しておくと良いでしょう。下の例では、上の４行をスキップ。列が、text（文字） か、numeric（数値）かを選択。読み込まないときには、skip しています。

df_un_fertility <- read_excel("./data/Age-specific fertility rates, Total fertility and .xls", 
    col_types = c("text", "numeric", "text", 
        "numeric", "text", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "skip", "text", "skip", 
        "skip", "skip", "skip", "skip", "skip"), 
    skip = 4)

df_un_fertility

Step 5. 列名を修正

修正したい列名をコピーしておく。

Country, ISO code, Period, Reference, NA, Total fertility

1列目から、6列目と、14列目と、15列目を、書き換えます。

colnames(df_un_fertility)[c(1:6,14:15)] <- c("country", "iso",  "period","year", "range",   "fertility_rate","mean_age","source")
df_un_fertility

Step 6. 日本のデータを確認

df_un_fertility_jp <- df_un_fertility %>% filter(country == "Japan")
df_un_fertility_jp

Step 7. 年齢ごとの数を比較しやすいように、データを変形

すこし、難しいですが、pivot_longer を Help で調べてください。

df_un_fertility_jp %>% pivot_longer(-c(1:6, 14:15), names_to = "age_range", values_to = "value")

Step 8. グラフを書いてみる。

df_un_fertility_jp %>% pivot_longer(-c(1:6, 14:15), names_to = "age_range", values_to = "value") %>%
  ggplot(aes(x=year, col=age_range, linetype = age_range)) + geom_line(aes(y=value)) 

いろいろなことがわかります。

Step 9. 出生平均年齢の表を降年齢順に見てみる。

df_un_fertility %>% 
  filter(year %in% c(2010)) %>% select(country, mean_age) %>% 
  drop_na(mean_age) %>%
  arrange(desc(mean_age))

Step 10. 密度分布グラフを書いてみる。

df_un_fertility %>% filter(year %in% c(1970, 1990, 2010)) %>% 
                                drop_na(mean_age) %>%
  ggplot(aes(mean_age, fill = factor(year))) + geom_density(alpha = 0.5) + 
#  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
#  coord_cartesian(xlim=c(20,40)) + 
  labs(title="Mean age at childbearing", fill = "year")

df_un_fertility %>% drop_na(mean_age) %>% group_by(year) %>% summarise(n=n()) %>% arrange(desc(n))

df_un_fertility %>% filter(year %in% c(1970, 1985, 1995, 2005, 2010)) %>% 
                                drop_na(mean_age) %>%
  ggplot(aes(mean_age, fill = factor(year))) + geom_density(alpha = 0.3) + 
  labs(title="Mean age at childbearing", fill = "year")

df_un_fertility %>% filter(year %in% c(1970, 1985,1995, 2005,2010)) %>% 
                                drop_na(mean_age) %>%
  ggplot(aes(mean_age, fill = factor(year))) + geom_density(alpha = 0.5) + facet_wrap(~ factor(year)) +
  labs(title="Mean age at childbearing", fill = "year")

11. 出生率の表を昇順に見てみる。

df_un_fertility %>% 
  filter(year %in% c(2010)) %>% select(country, fertility_rate) %>% 
  drop_na(fertility_rate) %>%
  arrange(fertility_rate)

Step 12. 密度分布グラフを書いてみる。

df_un_fertility %>% filter(year %in% c(1970, 1990, 2010)) %>% 
                                drop_na(mean_age) %>%
  ggplot(aes(fertility_rate, fill = factor(year))) + geom_density(alpha = 0.5) + 
  labs(title="Total fertility", fill = "year")

df_un_fertility %>% drop_na(fertility_rate) %>% group_by(year) %>% summarise(n=n()) %>% arrange(desc(n))

df_un_fertility %>% filter(year %in% c(1970, 1985, 1995, 2005, 2010)) %>% 
                                drop_na(mean_age) %>%
  ggplot(aes(fertility_rate, fill = factor(year))) + geom_density(alpha = 0.5) + 
  labs(title="Total fertility", fill = "year")

df_un_fertility %>% filter(year %in% c(1970, 1985, 1995, 2005, 2010)) %>% 
                                drop_na(mean_age) %>%
  ggplot(aes(fertility_rate, fill = factor(year))) + geom_density(alpha = 0.5) + facet_wrap(~ factor(year)) +
  labs(title="Total fertility", fill = "year")

3.2.1.1 参考例１：World Inequality Report 2022

WIR2022: https://ds-sl.github.io/data-analysis/wir2022.nb.html

3.2.1.2 参考例２:学校種類別進学率の推移(日本語データ)

男女共同参画局の資料です。

学校種類別進学率の推移: https://empowerment.tsuda.ac.jp/detail/82584

url_school_jp <- "https://www.gender.go.jp/about_danjo/whitepaper/r02/zentai/html/honpen/csv/zuhyo01-04-01.csv"

エンコーディング（Encoding type）を推測することができます。

guess_encoding(url_school_jp, n_max = 10000, threshold = 0.2)

df_school_jp <- read_csv(url_school_jp, locale = locale(encoding = "Shift_JIS"), skip=2)

Rows: 70 Columns: 10── Column specification ──────────────────────────────────────────────
Delimiter: ","
chr (1): 年度
dbl (9): 高等学校等（男子）, 高等学校等（女子）, 専修学校（専門課程，男子）, 専修学校（専門課程，女子）, 大学（学部...
ℹ 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_school_jp 

df_edu0 <- df_school_jp 
colnames(df_edu0) <- c("year", "highschool_m", "highschool_f", "vocational_m", "vocational_f", "university_m", "university_f", "juniorcollege_f", "gradschool_m", "gradschool_f")
df_edu00 <- df_edu0 %>% mutate(year = 1950:2019, 
                   highschool = (highschool_m + highschool_f)/2,
                   vocational = (vocational_m + vocational_f)/2,
                   university = (university_m + university_f)/2, 
                   juniorcollege = juniorcollege_f,
                   gradschool = (gradschool_m + gradschool_f)/2) 
df_edu00 %>% filter(year >= 1954) %>% select(-(2:10)) %>% 
  pivot_longer(3:5, names_to = "schools", values_to = "percentage") %>%
  mutate(types = factor(schools, levels = c("vocational", "juniorcollege", "university"))) %>%
  pivot_longer(c(highschool, gradschool), names_to = "highgrad", values_to ="value") %>%
  mutate(high_grad = factor(highgrad, levels = c("highschool", "gradschool"))) %>%
  ggplot() +
  geom_area(aes(x = year, y = percentage, fill = types)) +
  geom_line(aes(x = year, y = value, linetype = high_grad)) + 
  scale_x_continuous(breaks = round(seq(1960, 2020, by =10),1)) + 
  scale_y_continuous(breaks = round(seq(0, 100, by =10),1)) + 
  labs(title = "Tertially Education After Highschool", 
       subtitle = "with Highschool Graduates and Graduate School", 
       fill = "", linetype = "")

3.3 （続）探索的データ解析 (EDA)　

image from r4ds

NY.GDP.PCAP.CD: GDP per capita (current US$)

df_wdi_gdppcap <- WDI(country = "all", indicator = c(gdp_pcap = "NY.GDP.PCAP.CD"))
write_csv(df_wdi_gdppcap, "./data/df_wdi_gdppcap.csv")

df_wdi_gdppcap <- read_csv("./data/df_wdi_gdppcap.csv")

Rows: 16492 Columns: 5── Column specification ──────────────────────────────────────────────
Delimiter: ","
chr (3): country, iso2c, iso3c
dbl (2): year, gdp_pcap
ℹ 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_wdi_gdppcap

3.3.1 Data Transformation

3.3.1.1 列を select

df_wdi_gdppcap_small <- df_wdi_gdppcap %>% 
  select(country, year, gdp_pcap)
df_wdi_gdppcap_small

3.3.1.2 行を filter

df_wdi_gdppcap_short <- df_wdi_gdppcap %>% 
  filter(country %in% c("Japan", "Germany", "United States"))
df_wdi_gdppcap_short

df_wdi_gdppcap_small_short <- df_wdi_gdppcap %>% select(country, year, gdp_pcap) %>%
  filter(country %in% c("Japan", "Germany", "United States"))
df_wdi_gdppcap_small_short

次は、よく生じる、誤りの例で、ノコギリの歯（sawtoothed）のようなギザギザ・グラフと呼ばれます。なぜこのようなことが起きているかわかりますか。

df_wdi_gdppcap_small_short %>%
  ggplot(aes(x = year, y = gdp_pcap)) + geom_line()

df_wdi_gdppcap_small_short %>% filter(country %in% c("Japan")) %>%
  ggplot(aes(x = year, y = gdp_pcap)) + geom_line()

df_wdi_gdppcap_small_short %>%
  ggplot(aes(x = year, y = gdp_pcap)) + geom_point()

df_wdi_gdppcap_small_short %>% drop_na(gdp_pcap) %>%
  ggplot(aes(x = year, y = gdp_pcap, col = country)) + geom_line()

df_wdi_gdppcap_small_short %>% drop_na(gdp_pcap) %>%
  ggplot(aes(x = year, y = gdp_pcap, col = country)) + geom_line() +
  geom_point()

df_wdi_gdppcap_small_short %>% drop_na(gdp_pcap) %>%
  ggplot(aes(x = year, y = gdp_pcap)) + 
  geom_point(aes(color = country)) + 
  geom_smooth(method = 'loess', formula = 'y ~ x')

4 参考 - 今後の学習のために

4.1 RNotebook の活用

下のリンクを開き、右上の Code ボタンから、Download Rmd を選択すると、ダウンロードできますから、ダインロードしたものを、プロジェクト・フォールダーに移動またはコピーしてください。ダウンロードできないときは、Ctrl を押しながら、Download Rmd をクリックすると、Save As で保存できると思います。ブラウザーによって仕様が異なりますから、適切な方法を選んでください。

• https://ds-sl.github.io/intro2r/RNotebook-J.nb.html
• https://ds-sl.github.io/intro2r/Rmarkdown-J.nb.html

Windows でも、Mac でも提供されている、Google Chrome の場合には、Code ボタンから、ダンロードされるはずです。

4.2 クラウド - Posit Cloud

RStudio Cloudは、誰でもオンラインでデータサイエンスを行い、共有し、教え、学ぶことができる、軽量でクラウドベースのソリューションです。月25時間の制限がありますが、内容を共有して、他のアカウントから利用することも可能です。

クラウドサービス　How to Start Posit Cloud

1. Go to https://posit.cloud/
2. Sign Up: top right
3. Email address or Google account
4. New Project: Project Name

4.3 練習問題 Posit Primers

Posit Primers https://posit.cloud/learn/primers

最初の演習

Moodle（2月22日）にリンクをつけてあります。

• Visualization Basics
• Programming Basics

4.4 参考文献 References

• R For Data Science, by H. Wickham: https://r4ds.had.co.nz

  • Introduction: https://r4ds.had.co.nz/explore-intro.html#explore-intro

• Bookdown: https://bookdown.org, Archive

• Get Started: R Studio で R をはじめよう、R Markdown

• Introducton to R

• Data Analysis for Researchers 2022

• データ・サイエンスを始めましょう - Data Science for All