- ネタ元
- ノーベル賞受賞者データの取得
- スイスのチョコレート消費量データの入力
- ドイツ語による国名表記の英語への変換
- 国別チョコレート消費量データの取得
- 相関係数を算出
- 関係性をグラフで描画
- BMIデータを追加
- チョコレートの消費量とノーベル賞受賞者数、BMIの関係をグラフで描画
- チョコレート消費量とノーベル賞受賞者数の線形回帰モデル
■ネタ元
Chocolate and the Nobel Prize – a true story? » G-Forgeのコードを、ほぼそのまま借用しています。
■ノーベル賞受賞者データの取得
元記事のXMLパッケージ + readHTMLTable() 関数では動かなかったので、rvestパッケージを使ったコードに直しました。また、テーブルの列名などが少し変更されていたので対応しました。
library(rvest) # 元記事のXML + readHTMLでは動かなかったので theurl <- "https://en.wikipedia.org/wiki/List_of_countries_by_Nobel_laureates_per_capita" tables <- read_html(theurl) %>% html_nodes("table") %>% .[1] %>% html_table(trim = T) nobel_prizes <- tables[[1]] # Clean column names colnames(nobel_prizes) <- gsub(" ", "_", gsub("(/|\\[[0-9]+\\])", "", gsub("\n", " ", colnames(nobel_prizes)) ) ) # Delete those that aren't countries and thus lack rank nobel_prizes$Rank <- as.numeric(as.character(nobel_prizes$Rank)) nobel_prizes <- subset(nobel_prizes, is.na(Rank) == FALSE) # Clean the country names nobel_prizes$Entity <- gsub("[^a-zA-Z ]", "", nobel_prizes$Entity) # Clean the loriates variable nobel_prizes$Laureates10_million <- as.numeric(as.character(nobel_prizes$Laureates10_million))
■スイスのチョコレート消費量データの入力
スクレイピング先にデータがないため、手入力する必要があるようです。ここも、列名 (Country → Entity) の変更があります。
nobel_prizes$Chocolate_consumption <- NA # http://www.chocosuisse.ch/web/chocosuisse/en/documentation/faq.html nobel_prizes$Chocolate_consumption[nobel_prizes$Entity == "Switzerland"] <- 11.9
■ドイツ語による国名表記の英語への変換
スクレイピング先のサイトがドイツ語で書かれているためです。rvestを用いた記述に書き換えました。
# Translation from German to English theurl <- "https://www.thoughtco.com/countries-of-the-world-index-4101906" tables <- read_html(theurl) %>% html_nodes("table") %>% .[1] %>% html_table(trim = T) translate_german <- tables[[1]] translate_german <- translate_german[3:nrow(translate_german), 1:2] colnames(translate_german) <- c("English", "German") translate_german$German <- gsub("([ ]+(f|pl)\\.$|\\([[:alnum:] -]+\\))", "", translate_german$German)
■国別チョコレート消費量データの取得
ドイツのチョコレートメーカー? のサイトに記載されている、国別のチョコレート消費量を取得します。rvestを用いた記述に書き換えました。また、ノーベル賞受賞者数のテーブルではCountryではなく、Entity (地域?) になっているため、最後のSQL文を書き換えました。
# Get the consumption from a German list theurl <- "http://www.theobroma-cacao.de/wissen/wirtschaft/international/konsum" tables <- read_html(theurl) %>% html_nodes("table") %>% .[1] %>% html_table(trim = T) german_chocolate_data <- tables[[1]][2:nrow(tables[[1]]), ] names(german_chocolate_data) <- c("Country", "Chocolate_consumption") german_chocolate_data$Country <- gsub("([ ]+f\\.$|\\([[:alnum:] -]+\\))", "", german_chocolate_data$Country) library(sqldf) sql <- paste0("SELECT gc.*, tg.English as Country_en", " FROM german_chocolate_data AS gc", " LEFT JOIN translate_german AS tg", " ON gc.Country = tg.German", " OR gc.Country = tg.English") german_chocolate_data <- sqldf(sql) german_chocolate_data$Country <- ifelse(is.na(german_chocolate_data$Country_en), german_chocolate_data$Country, german_chocolate_data$Country_en) german_chocolate_data$Country_en <- NULL german_chocolate_data$Chocolate_consumption_tr <- NA for (i in 1:nrow(german_chocolate_data)) { number <- as.character(german_chocolate_data$Chocolate_consumption[i]) if (length(number) > 0) { m <- regexpr("^([0-9]+,[0-9]+)", number) if (m > 0) { german_chocolate_data$Chocolate_consumption_tr[i] <- as.numeric( sub(",", ".", regmatches(number, m)) ) } else { m <- regexpr("\\(([0-9]+,[0-9]+)", number) if (m > 0) german_chocolate_data$Chocolate_consumption_tr[i] <- as.numeric( sub("\\(", "", sub(",", ".", regmatches(number, m)) ) ) } } } sql <- paste0("SELECT np.*, gp.Chocolate_consumption_tr AS choc", " FROM nobel_prizes AS np", " LEFT JOIN german_chocolate_data AS gp", " ON gp.Country = np.Entity") nobel_prizes <- sqldf(sql)
■相関係数を算出
res_cor <- nobel_prizes %>% filter(!is.na(choc)) %>% select(Laureates10_million, choc) %>% cor() res_cor
Laureates10_million choc Laureates10_million 1.0000000 0.7336186 choc 0.7336186 1.0000000
■関係性をグラフで描画
一部、列名が現在取得できるテーブルのものと異なっていたので修正しました。
library(ggplot2) ggplot(data = subset(nobel_prizes, is.na(choc) == FALSE), aes(x = choc, y = Laureates10_million)) + ylab("1千万人あたりのノーベル賞受賞者数") + xlab("1人あたりのチョコレート消費量 (kg)") + geom_point(size = 4, colour = "#444499") + geom_text(aes(label = Entity), size = 5, position = position_jitter(width = 0.5, height = 0.5), colour = "#990000") + annotate("text", x = 1.8, y = 25, label=paste("相関係数:", round(res_cor[1, 2], 2)), size = 6)

■BMIデータを追加
ついでに、BMI (Body mass index) の国別データも取得し、追加してみたそうです。OECDのライブラリから、2012年の調査データを取得するようですが、URLが変更されていたので、修正しました。
# Percentage of males with a BMI > 25 kg/m2 tables <- read_html("https://www.oecd-ilibrary.org/sites/ovob-ml-table-2012-2-en/index.html?itemId=/content/component/ovob-ml-table-2012-2-en") %>% html_nodes("table") %>% .[1] %>% html_table(trim = T) ob <- tables[[1]] ob[, 2] <- as.character(ob[, 2]) ob <- apply(ob, FUN = as.character, MARGIN = 2) ob <- ob[, 2:ncol(ob)] ob <- ob[4:nrow(ob), ] last_obesitas <- apply(apply(ob[, 2:ncol(ob)], FUN = as.numeric, MARGIN = 2), MARGIN = 1, FUN = function(x) { if (any(is.na(x) == FALSE)) return(x[max(which(is.na(x) == FALSE))]) else return(NA) }) ob <- data.frame(Country = ob[, 1], last_obesitas = last_obesitas) ob <- subset(ob, is.na(last_obesitas) == FALSE) sql <- paste0("SELECT np.*, ob.last_obesitas AS obesitas", " FROM nobel_prizes AS np", " LEFT JOIN ob", " ON ob.Country = np.Entity") nobel_prizes <- sqldf(sql)
■チョコレートの消費量とノーベル賞受賞者数、BMIの関係をグラフで描画
ggplot(data = subset(nobel_prizes, is.na(choc) == FALSE), aes(x = choc, y = Laureates10_million)) + ylab("1千万人あたりのノーベル賞受賞者数") + xlab("1人あたりのチョコレート消費量 (kg)") + geom_point(aes(size = obesitas), colour = "#444499") + scale_size(range = c(2, 10)) + geom_text(aes(label = Entity), size = 5, position = position_jitter(width = 0.5, height = 0.5), colour = "#770000")

■チョコレート消費量とノーベル賞受賞者数の線形回帰モデル
目的変数を受賞者数 (の割合)、説明変数をチョコレートの消費量としているので、「ノーベル賞の受賞者数は、チョコレートの消費量によって説明できる」と仮定した因果モデルですね。列名が少し異なっているので修正しました。
fit <- lm(Laureates10_million ~ choc, data = nobel_prizes) summary(fit)
Call: lm(formula = Laureates10_million ~ choc, data = nobel_prizes) Residuals: Min 1Q Median 3Q Max -12.552 -3.884 -1.087 2.534 20.362 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -1.1304 3.5830 -0.316 0.75646 choc 2.4593 0.5695 4.318 0.00053 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 7.429 on 16 degrees of freedom (56 observations deleted due to missingness) Multiple R-squared: 0.5382, Adjusted R-squared: 0.5093 F-statistic: 18.65 on 1 and 16 DF, p-value: 0.0005302