まず必要なパッケージのダウンロードやデータの読み込みを行った。データの読み込みについては、全体を通してread_csv関数を用いて行った。
折れ線グラフについては、qplotを用いて、年を色分けして作成した。その際colour=yearであると、年の違いがグラデーションで表され、線もほぼ一本になってしまったことから、colour=factor(year)によって数値を文字化し色を分け、group=yearによって年でグループ分けして線が一つになってしまうのを防いだ。
library("tidyverse")
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.3 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library("ggplot2")
library("magrittr")
##
## Attaching package: 'magrittr'
##
## The following object is masked from 'package:purrr':
##
## set_names
##
## The following object is masked from 'package:tidyr':
##
## extract
heatstroke <- read_csv("heatstroke_2021_25.csv")
## Rows: 25 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (3): year, month, heatstroke_emergency_transports
## date (1): date
##
## ℹ 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.
heatstroke %<>% rename(月別熱中症救急搬送者件数='heatstroke_emergency_transports')
qplot(month, 月別熱中症救急搬送者件数, colour=factor(year),group=year,data=heatstroke,geom="line")
## Warning: `qplot()` was deprecated in ggplot2 3.4.0.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
スクレイピングについては、授業内の方法と同様にXPATHを取得して行った。収集したデータは縦型でなかったため、1列目をyearと名付け、yearと必要な月を取り出したデータを作成した。次にpivot_longer関数を使ってデータを縦型にして、数値化などの作業を行った。
グラフについては、一定期間での支出の変化を観察することが目的のグラフであることから、(1)同様折れ線グラフで表した。作業内容に関しても、(1)と違いはない。
library("rvest")
##
## Attaching package: 'rvest'
## The following object is masked from 'package:readr':
##
## guess_encoding
url <- "https://www.icecream.or.jp/iceworld/data/expenditures.html"
ice <- read_html(url) %>%
html_element(xpath="/html/body/div[1]/div/div[3]/div/div[2]/div/div/div[4]/table") %>% html_table
library(dplyr)
colnames(ice)[1] <- "year"
ice %>% select(year,"5月","6月","7月","8月","9月") -> ice1
ice1 %<>% pivot_longer(cols=c("5月", "6月", "7月", "8月", "9月"),names_to="month",values_to="アイスクリーム月別支出金額")
ice1 %<>% mutate(アイスクリーム月別支出金額=parse_number(アイスクリーム月別支出金額))
qplot(month, アイスクリーム月別支出金額, colour=factor(year),group=year,data=ice1,geom="line")
散布図を作成するにあたり、一つのデータにまとめる必要があったため、左外部結合を行った。heatstrokeデータに関しては、year、monthの列に単位ついていなかったことから、paste0関数を用いて年、月の単位をつけ、結合ができる形にした。
散布図については、ggplotを用いて作成し、色を年で分け、点が月で表示されるようにgeom_point形式からgeom_text形式に変更した。また支出金額のすべての点が示されるように、横軸の表示範囲を調整した。
散布図からは、アイスクリーム月別支出金額と月別熱中症救急搬送者件数の間に強い正の相関がみられる。内訳については、最も気温が高いと推測される7月、8月が図右上に多いことがわかり、気温が比較的低い月は右下に集まっている。
しかし両者には直接的な因果関係はないので、気温の高さが熱中症の増加に影響を及ぼし、またアイスの売り上げにも影響を与えるという、気温を介在とした関係があると思われる。
heatstroke %>% select(year, month, 月別熱中症救急搬送者件数) ->heatstroke1
heatstroke1$year <- paste0(heatstroke1$year,"年")
heatstroke1$month <- paste0(heatstroke1$month, "月")
ice1 %>% left_join(heatstroke1,by=c("year","month")) ->sanpuzu
ggplot(data=sanpuzu, mapping=aes(x=`アイスクリーム月別支出金額`,y=`月別熱中症救急搬送者件数`,label=month, colour=year)) + geom_text() + scale_x_continuous(limits=c(800,2000))
まず散布図を描いて、lm関数を用いて回帰分析を行った。次に回帰直線や決定係数の出力、回帰分析の診断グラフの表示を行った。
決定係数は0.85を超えているため、あてはまりが強く、アイスへの支出が熱中症搬送者の人数に強い関係を持つと考えられる。しかし診断グラフの一つ目から、赤線が水平にあまり近くないことから、直線の当てはまりがあまり妥当でないことが考えられる。二つ目のグラフからはモデルのあてはめが、三つ目のグラフからモデルの前提が正しいことがうかがえる。
総じて説明変数と目的変数の間には強い関係があると考えられるが、直接の因果関係が存在しないことには留意する必要がある。
sanpuzu %$% plot(アイスクリーム月別支出金額,月別熱中症救急搬送者件数)
lm(月別熱中症救急搬送者件数~アイスクリーム月別支出金額,sanpuzu) -> kaiki
kaiki
##
## Call:
## lm(formula = 月別熱中症救急搬送者件数 ~ アイスクリーム月別支出金額,
## data = sanpuzu)
##
## Coefficients:
## (Intercept) アイスクリーム月別支出金額
## -40353.35 42.06
kaiki %>% abline(col="red")
kaiki %>% summary()
##
## Call:
## lm(formula = 月別熱中症救急搬送者件数 ~ アイスクリーム月別支出金額,
## data = sanpuzu)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8291.7 -3116.4 -507.7 2323.5 12715.2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -40353.350 4882.708 -8.265 2.44e-08 ***
## アイスクリーム月別支出金額 42.062 3.541 11.877 2.71e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5143 on 23 degrees of freedom
## Multiple R-squared: 0.8598, Adjusted R-squared: 0.8537
## F-statistic: 141.1 on 1 and 23 DF, p-value: 2.713e-11
kaiki %>% plot()
csvの読み込みについて、ファイルに文字化けが起きていたり、複数行にわたって列名が入っていたため、列名の割り当てや文字コードの修正を行った(この点については、生成AIを参考にしながら作業を行った)。
次に列名や単位を変更して、すでに作ったデータと結合を行った。回帰分析に方法については、(4)と同様に進めた。
決定係数の当てはまりは0.6台でまずまずの結果を示している。一方、回帰直線の係数は88と(4)の回帰直線より高い。診断グラフの結果から、モデルの当てはまりなどは高いことが分かった。
以上のことから、最高気温はアイスへの支出に大きな影響を与えているといえる。決定係数の高さから一定の説明力が担保されており、今までの疑似相関と異なり、因果関係が認められる。
kisyoumoto <- read_csv("weather_sum_2021-25.csv")
## Warning: One or more parsing issues, call `problems()` on your data frame for details,
## e.g.:
## dat <- vroom(...)
## problems(dat)
## Rows: 29 Columns: 1
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): ダウンロードした時刻:2026/07/15 16:26:33
##
## ℹ 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.
kisyou <- read_csv("weather_sum_2021-25.csv", skip = 6,
col_names = FALSE,
locale = locale(encoding = "UTF-8"))
## Rows: 25 Columns: 24
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): X1
## dbl (23): X2, X3, X4, X5, X6, X7, X8, X9, X10, X11, X12, X13, X14, X15, X16,...
##
## ℹ 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.
kisyou1 <- kisyou %>% select(X1, X16) %>%
rename(YM = X1, 月ごとの平均最高気温 = X16) %>%
separate(YM, into = c("year", "month"), sep = "/") %>%
mutate(year = paste0(year,"年"),month = paste0(month,"月"))
sanpuzu %>% left_join(kisyou1,by=c("year","month")) -> sanpuzu2
sanpuzu2 %$% plot(月ごとの平均最高気温,アイスクリーム月別支出金額)
lm(アイスクリーム月別支出金額~月ごとの平均最高気温,sanpuzu2) ->kaiki2
kaiki2
##
## Call:
## lm(formula = アイスクリーム月別支出金額 ~ 月ごとの平均最高気温,
## data = sanpuzu2)
##
## Coefficients:
## (Intercept) 月ごとの平均最高気温
## -1686.56 88.11
kaiki2 %>% abline(col="blue")
kaiki2 %>%summary()
##
## Call:
## lm(formula = アイスクリーム月別支出金額 ~ 月ごとの平均最高気温,
## data = sanpuzu2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -411.49 -123.95 9.04 84.22 333.51
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1686.56 446.25 -3.779 0.000971 ***
## 月ごとの平均最高気温 88.11 12.92 6.820 5.9e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 174.2 on 23 degrees of freedom
## Multiple R-squared: 0.6692, Adjusted R-squared: 0.6548
## F-statistic: 46.52 on 1 and 23 DF, p-value: 5.903e-07
kaiki2 %>% plot()
(ⅰ)と同様の作業を行った。
決定係数は0.6台で、係数は3000台であるが、縦軸の単位から(ⅰ)と同様の影響力とみることができる。診断グラフから見たデータやモデルの当てはまりは強い。
よって最高気温は熱中症による搬送者数へも強い影響を与えており、因果関係が存在すると考えられる、
sanpuzu2 %$% plot(月ごとの平均最高気温,月別熱中症救急搬送者件数)
lm(月別熱中症救急搬送者件数~月ごとの平均最高気温,sanpuzu2) ->kaiki3
kaiki3
##
## Call:
## lm(formula = 月別熱中症救急搬送者件数 ~ 月ごとの平均最高気温,
## data = sanpuzu2)
##
## Coefficients:
## (Intercept) 月ごとの平均最高気温
## -115925 3840
kaiki3 %>% abline(col="green")
kaiki3 %>%summary()
##
## Call:
## lm(formula = 月別熱中症救急搬送者件数 ~ 月ごとの平均最高気温,
## data = sanpuzu2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16404 -5069 -492 6563 15873
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -115925.2 21755.6 -5.329 2.07e-05 ***
## 月ごとの平均最高気温 3840.4 629.8 6.098 3.21e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8491 on 23 degrees of freedom
## Multiple R-squared: 0.6179, Adjusted R-squared: 0.6012
## F-statistic: 37.19 on 1 and 23 DF, p-value: 3.208e-06
kaiki3 %>% plot()
気温と温度から求められる「蒸し暑さ」の指数である、「不快指数」を新たに変数として用意して、アイスクリームへの支出と熱中症患者数との関係を調べる。ちなみに日本人は不快指数85で不快感を感じるとされているが、風向きや風速などで条件が変化するため、気象庁に公式で使われている指数ではない。
不快指数は、0.81×気温+0.01×湿度×(0.99×気温ー14.3)+46.3 という式で求められる。
気象についてのデータから、平均気温と平均湿度のデータを利用して、mutate関数を用いて不快指数を導く変数を新たに作成した。結合ができるように列名や単位を変更し、同様の方法で回帰分析を行った。
不快指数とアイスクリームの売り上げの関係については、決定係数が0.7を超えていて、あてはまりがかなり強いことが分かった。回帰直線の傾きから、不快指数が支出に一定程度の影響を与えていることが読み取れる。
不快指数と熱中症による搬送者数の関係は、決定係数は0.6であるが、回帰直線の傾きは正であることから、不快指数が熱中症の発生にもある程度影響を与えていることが分かる。
どちらのグラフからも、不快指数が上がれば支出、搬送者数が上がるという関係を導くことができる。不快指数は気温や湿度を媒介にした指数であることから、気温だけを説明変数とするよりも正しい結果を導くと考えられる。今回のような複数のデータから作り出された指数によって予測精度が向上したならば、より多くのデータを盛り込んだ指数があれば、より正確な結果を導くことができると推測される。
kisyou %<>% mutate(X23 = 0.81*X2+0.01*X13*(0.99*X2-14.3)+46.3)
hukai <- kisyou %>% select(X1, X23) %>%
rename(YM = X1, 不快指数 = X23) %>%
separate(YM, into = c("year", "month"), sep = "/") %>%
mutate(year = paste0(year,"年"),month = paste0(month,"月"))
sanpuzu %>% left_join(hukai,by=c("year","month")) -> sanpuzu3
sanpuzu3 %$% plot(不快指数,アイスクリーム月別支出金額)
lm(アイスクリーム月別支出金額~不快指数,sanpuzu3) ->kaiki4
kaiki4
##
## Call:
## lm(formula = アイスクリーム月別支出金額 ~ 不快指数,
## data = sanpuzu3)
##
## Coefficients:
## (Intercept) 不快指数
## -2034.08 45.41
kaiki4 %>% abline(col="blue")
kaiki4 %>%summary()
##
## Call:
## lm(formula = アイスクリーム月別支出金額 ~ 不快指数,
## data = sanpuzu3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -257.00 -101.32 36.74 101.24 258.16
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2034.08 402.57 -5.053 4.09e-05 ***
## 不快指数 45.41 5.39 8.424 1.75e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 149.8 on 23 degrees of freedom
## Multiple R-squared: 0.7552, Adjusted R-squared: 0.7446
## F-statistic: 70.96 on 1 and 23 DF, p-value: 1.748e-08
kaiki4 %>% plot()
sanpuzu3 %$% plot(不快指数,月別熱中症救急搬送者件数)
lm(月別熱中症救急搬送者件数~不快指数,sanpuzu3) ->kaiki5
kaiki5
##
## Call:
## lm(formula = 月別熱中症救急搬送者件数 ~ 不快指数,
## data = sanpuzu3)
##
## Coefficients:
## (Intercept) 不快指数
## -123326 1875
kaiki5 %>% abline(col="blue")
kaiki5 %>%summary()
##
## Call:
## lm(formula = 月別熱中症救急搬送者件数 ~ 不快指数,
## data = sanpuzu3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18597 -4831 1373 4552 15464
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -123326.2 22574.0 -5.463 1.49e-05 ***
## 不快指数 1875.2 302.2 6.204 2.49e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8401 on 23 degrees of freedom
## Multiple R-squared: 0.626, Adjusted R-squared: 0.6097
## F-statistic: 38.49 on 1 and 23 DF, p-value: 2.493e-06
kaiki5 %>% plot()
この授業を通して、データを読み解く能力、そして活用する能力が身についたと思う。特に後半のRに関しては、実践的なデータサイエンスの力を養うことができた。情報化社会を生き抜くために、学んだことを生かしながら、もっと多くの知識をつけていきたい。