大学入試センター令和7年度試験の問題作成の方向性,試作問題等のなかで、情報Ⅰデータサイエンス分野の参考問題が掲載されています。
これらの要請にかなう問題として作られ、データは公開されている実データが使われています。
日頃からRやPythonなどを利用したデータ分析に親しんでおくことが求められています。
データは、
一般社団法人日本冷凍空調工業会
一般社団法人日本アイスクリーム協会
気象庁
(金沢市のデータ)
からダウンロードまたはコピーし、エクセルで整形しました。
分析に必要なパッケージを読み込みます。
library(TSstudio)
library(MLmetrics)
##
## 次のパッケージを付け加えます: 'MLmetrics'
## 以下のオブジェクトは 'package:base' からマスクされています:
##
## Recall
library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(TTR)
library(quantmod)
## 要求されたパッケージ xts をロード中です
## 要求されたパッケージ zoo をロード中です
##
## 次のパッケージを付け加えます: 'zoo'
## 以下のオブジェクトは 'package:base' からマスクされています:
##
## as.Date, as.Date.numeric
##
## 次のパッケージを付け加えます: 'xts'
## 以下のオブジェクトは 'package:dplyr' からマスクされています:
##
## first, last
エアコンの売上台数とアイスクリームの売上個数のデータを読み込みます。
気温(kion)、湿度(situdo)は問5で使用します。
x <- read.csv("data/ice_aircon.csv")
str(x)
## 'data.frame': 60 obs. of 5 variables:
## $ date : chr "2016_1" "2016_2" "2016_3" "2016_4" ...
## $ aircon: int 434 504 769 420 759 1470 1542 651 469 336 ...
## $ ice : int 464 397 493 617 890 883 1292 1387 843 621 ...
## $ kion : num 5 4.9 8.3 14.1 19.4 22.4 26.2 27.5 23.8 18 ...
## $ situdo: int 74 71 63 64 64 71 74 71 78 68 ...
head(x)
date | aircon | ice | kion | situdo |
---|---|---|---|---|
2016_1 | 434 | 464 | 5.0 | 74 |
2016_2 | 504 | 397 | 4.9 | 71 |
2016_3 | 769 | 493 | 8.3 | 63 |
2016_4 | 420 | 617 | 14.1 | 64 |
2016_5 | 759 | 890 | 19.4 | 64 |
2016_6 | 1470 | 883 | 22.4 | 71 |
売上の変化を図示します。
par(mfrow = c(2, 1))
plot(x$aircon, type='l',ylim=c(0,1700))
plot(x$ice, type='l' ,ylim=c(0,1700))
airconとiceのグラフを見ると形は似通っていますが、iceのピークが少し右にずれているようです。
airconとiceは気温の影響を受けます。この場合、airconとiceには気温による擬似相関が発生しますが、このことは問5に出てきます。
iceの売上は、その日の気温に即時的に反応し、airconは前もって準備するので
airconのピークが少し前に来ると考えられます。
分析のため、データを時系列データに変換します。
aircon<-ts(x[,2],start=c(2016,1),frequency=12)
ice<-ts(x[,3],start=c(2016,1),frequency=12)
aircon
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2016 434 504 769 420 759 1470 1542 651 469 336 427 571
## 2017 520 547 816 482 882 1560 1591 714 495 317 435 565
## 2018 515 582 916 534 901 1568 1763 818 557 376 493 627
## 2019 570 608 1000 651 1078 1662 1227 960 731 324 419 585
## 2020 481 551 905 598 947 1844 1464 940 632 379 494 635
ice
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2016 464 397 493 617 890 883 1292 1387 843 621 459 562
## 2017 489 382 472 624 915 914 1394 1370 826 599 489 573
## 2018 507 416 607 746 894 1021 1506 1443 861 640 492 537
## 2019 494 423 542 667 1000 991 1236 1513 996 724 531 584
## 2020 510 482 610 689 1040 1123 1155 1658 1025 649 573 599
時系列データを図示すると、横軸は Time になります。
par(mfrow = c(2, 1))
ts.plot(aircon)
ts.plot(ice)
一般に、時系列データは、次の成分に分解されます。
観測値 (Obserbed) | 実際の値 | |
トレンド (Trend) | データの長期的変化・傾向 | |
季節変動 (Seasonal) | 周期変動 | |
誤差変動 (Random) | 誤差による変化 | |
トレンドは12ヶ月ごとの移動平均で算出されます。
関数を利用して時系列分解(time series decomposition)を行います。
ts_decompose(ice)
ts_decompose(aircon)
トレンドを見ると、
アイスクリームは順調に売上を伸ばしているようですが、エアコンは2020にやや落ち込んでいるように見えます。
誤差変動は2019末に大きくマイナスに振れており、何か異常事態が起きていることが察せられます。ちょうどコロナの始まりと重なっているようです。
移動平均を具体的に見てみましょう。
トレンドは12ヶ月移動平均なので、SMA
関数を用いて12ヶ月移動平均に変換します。
i.ice<-SMA(ice,12)
i.ice
## Jan Feb Mar Apr May Jun Jul Aug
## 2016 NA NA NA NA NA NA NA NA
## 2017 744.4167 743.1667 741.4167 742.0000 744.0833 746.6667 755.1667 753.7500
## 2018 755.4167 758.2500 769.5000 779.6667 777.9167 786.8333 796.1667 802.2500
## 2019 804.7500 805.3333 799.9167 793.3333 802.1667 799.6667 777.1667 783.0000
## 2020 809.7500 814.6667 820.3333 822.1667 825.5000 836.5000 829.7500 841.8333
## Sep Oct Nov Dec
## 2016 NA NA NA 742.3333
## 2017 752.3333 750.5000 753.0000 753.9167
## 2018 805.1667 808.5833 808.8333 805.8333
## 2019 794.2500 801.2500 804.5000 808.4167
## 2020 844.2500 838.0000 841.5000 842.7500
i.aircon<-SMA(aircon,12)
i.aircon
## Jan Feb Mar Apr May Jun Jul Aug
## 2016 NA NA NA NA NA NA NA NA
## 2017 703.1667 706.7500 710.6667 715.8333 726.0833 733.5833 737.6667 742.9167
## 2018 743.2500 746.1667 754.5000 758.8333 760.4167 761.0833 775.4167 784.0833
## 2019 808.7500 810.9167 817.9167 827.6667 842.4167 850.2500 805.5833 817.4167
## 2020 810.5000 805.7500 797.8333 793.4167 782.5000 797.6667 817.4167 815.7500
## Sep Oct Nov Dec
## 2016 NA NA NA 696.0000
## 2017 745.0833 743.5000 744.1667 743.6667
## 2018 789.2500 794.1667 799.0000 804.1667
## 2019 831.9167 827.5833 821.4167 817.9167
## 2020 807.5000 812.0833 818.3333 822.5000
移動平均なので、2016の12月になって初めて値が現れます。
par(mfrow = c(2, 1))
ts.plot(i.aircon)
ts.plot(i.ice)
12ヶ月移動平均のグラフはトレンドと一致しています。
移動平均がトレンドであることを知らないと、イメージしにくい問題です。
移動期間をいろいろ変えて図示してみます。
aircon
aircon9 <-SMA(aircon,9)
aircon6 <-SMA(aircon,6)
aircon3 <-SMA(aircon,3)
par(mfrow = c(3, 1))
ts.plot(aircon9)
ts.plot(aircon6)
ts.plot(aircon3)
⓪は3ヶ月移動平均 、②は9ヶ月移動平均 、③は6ヶ月移動平均
①は12ヶ月移動平均でトレンドになります。
自己相関係数
未来は過去の影響を受けます。
時系列のある時刻のデータとkだけ離れた時刻のデータとの相関係数を、ラグkの自己相関係数と言い、\(\rho_k\) と書きます。
(数理的な説明は別の機会に譲ります)
自己相関(ACF)と偏自己相関(PACF)
ts_cor(aircon, lag.max = 60)
自己相関関数(Autocorrelation
Function)は、1年ごとにピーク(赤)を迎え、時間の経過とともに小さくなっていきます。
偏自己相関(Partial AutoCorrelation
Function)は、ラグhに関して時点t-hと時点tの間に存在するh-1個の観測値の影響を除去したあとの相関ですが、難しいのでここでは立ち入りません。
問3のグラフは自己相関(ACF)のグラフです。
ラグごとの散布図
ts_lags(aircon, lags = 1:12)
ラグ1で正の相関であったものが、ラグ5では負の相関になり、ラグ12(1年後)では強い正の相関となっており、これがACFのグラフのはじめの赤い棒に相当します。
この変化は、12ヶ月ごとに繰り返されます。
相互相関(cross-correlation function)CCF
過去の自分との類似性を見るのが「自己相関」、他の時系列データとの類似性を見るのが「相互相関」です。
セールスアナリティクス様「時系列データの自己相関と相互相関をPythonで求めてみよう」参照
ccf(ice, aircon)
#相関係数
cor(ice, aircon)
## [1] 0.61977
Lag=0:ACF=0.62 は相関係数
Lag=1:ACF=0.86
は相互相関の最大値で、airconの売上変動を1ヶ月ずらすと、iceの売上変動に近いことがわかります。
重回帰分析
これまではairconの売上とiceの売上の関係を時系列で分析してきました。
しかし実際には、他にも複数の要因がありそれらの要因が影響しあっています。
複数の要因(変数)の関係を分析するのが多変量解析です。
airconとiceのほかにkion(気温)、situdo(湿度)を含めた4変数の分析を行います。
今回は時系列分析ではないので、1列目のdateを削除します。
x <- x[,-1]
str(x)
## 'data.frame': 60 obs. of 4 variables:
## $ aircon: int 434 504 769 420 759 1470 1542 651 469 336 ...
## $ ice : int 464 397 493 617 890 883 1292 1387 843 621 ...
## $ kion : num 5 4.9 8.3 14.1 19.4 22.4 26.2 27.5 23.8 18 ...
## $ situdo: int 74 71 63 64 64 71 74 71 78 68 ...
パッケージ “psych”(心理学)の関数 pairs.panelを用いて、散布図・相関行列を図示します。
library(psych)
##
## 次のパッケージを付け加えます: 'psych'
## 以下のオブジェクトは 'package:MLmetrics' からマスクされています:
##
## AUC
psych::pairs.panels(x)
相関が強いのは ice と kion であり、相関係数は0.91
と高くなっています。
擬似相関
散布図を描くと、airconとiceの間には相関があるように見えますが、気温が上がると、airconとiceともに増えることからきた擬似相関と言います。
airconの売上を目的変数として、ice,kion,situdo の3変数で重回帰分析を行ってみます。
(res<-lm(aircon~.,data=x))
##
## Call:
## lm(formula = aircon ~ ., data = x)
##
## Coefficients:
## (Intercept) ice kion situdo
## 393.167 1.060 -13.712 -3.508
summary(res)
##
## Call:
## lm(formula = aircon ~ ., data = x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -586.42 -184.16 -67.42 96.42 823.98
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 393.1669 607.0775 0.648 0.519865
## ice 1.0601 0.2977 3.561 0.000762 ***
## kion -13.7116 12.3127 -1.114 0.270201
## situdo -3.5076 8.9389 -0.392 0.696256
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 321.9 on 56 degrees of freedom
## Multiple R-squared: 0.4, Adjusted R-squared: 0.3678
## F-statistic: 12.44 on 3 and 56 DF, p-value: 2.41e-06
実際に影響力を持つ変数は”kion”のみです。
(アスタリスクのついたものが影響力を持つ変数です)
モデルの適合度を示す “Adjusted R-squared: 0.3678”
であり、3つの変数を使うのは良いモデルとは言えないようです。
回帰診断図 (diagnostic plots)
par(mfrow=c(2,2),oma = c(1,1,2,1),mar = c(4, 4, 2, 1))
plot(res,pch=21,bg=2,col=2,cex=1.5)
数字のついたデータは外れ値で、あまりモデルの当てはまりはよくないようです。
重回帰分析についてもっと詳しく知りたい方は、当サイト
問題は、将来のデータサイエンスの学びを念頭において作られています。
また、丁寧な誘導・説明があり、よく考えれば解答できるようになっています。
内容は、時系列分析、回帰分析など情報Ⅰの教科書にあるものですが、文章だけでは実感が湧きません。
統計ソフトを用いて実際に分析を体験しておくことが有用であると感じられます。