トップ 履歴 一覧 Farm ソース 検索 ヘルプ RSS ログイン

チョコレートの消費量とノーベル賞受賞者の関係


キーワード

最終更新時間：2020年06月02日 18時27分31秒
アフィリエイト・広告について
プライバシーポリシー

RPubsにも同じ記事があります。

■ネタ元

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

カテゴリ: [R,データ分析]

https://mana.bi/ltl_logo.png

Menu

最近更新したページ

RcppMeCabパッケージ / Jupyter Notebookでボタン1つで日本語PDFを作成する / ggimageパッケージの例_東京芝2400mのレコード推移 / データサイエンスメモ / チョコレートの消費量とノーベル賞受賞者の関係

アクセスの多いページ


Copyright (C) 2020 タナカケンタ
Powered by FreeStyleWiki3.6.5 with Perl5.026001