時系列予測モデルの評価指標について
以下の資料を参考に、時系列予測モデルの代表的な評価指標について、特徴や、メリットデメリットをまとめてみます。本記事の内容は、参考資料に記載されている内容から私がざっくり解釈したもののため、正確性には欠けていると思います。
https://otexts.com/fpp3/accuracy.html
まず、評価指標は3つの種類に分類されます。
- Scale-dependent Errors
- Percentage Errors
- Scaled Errors
Scale-dependent Errorsは予測値と実測値の差分から計算されるのもので、単位に依存します。そのため、異なる時系列データの予測モデルのScale-dependent Errors（例えば、株価を予想したAモデルのScale-dependent Errorsと金利を予想したBモデルのScale-dependent Errors）は比較できません。
この点を克服したのがPercentage Errorsです。実測値等を使用して、予測誤差を%で表します。
しかし、Percentage Errorsにも問題があり、誤差を％で表すために、予測誤差を実測値等で割る際、実測値が0の場合は定義できません。また、0でなくとも0に近い値で割ろうとすると、急に大きな値になってしまうため、安定しません。
Scaled ErrorsはPercentage ErrorsとScale-dependent Errorsの良いとこ取りを目指した指標です。
以下では具体的な指標例とともに、上記の内容を掘り下げます。
Scale-dependent Errors
まず予測誤差は以下のように表せます。
が実測値、
が予測値です。
ポピュラーなScale-dependent Errorsの指標として以下が挙げられます。
- Mean Absolute Error : MAE =
- Root Mean Squared Error : RMSE =
MAEは予測誤差の絶対値の平均を取ったものです。これを評価指標として用いてモデル選択を行うと、その予測モデルは中央値を予測するモデルになります。これはどういうことかというと、「MAEを評価指標として使う」=「MAEが最小となるモデルを選択する」ということになるのですが、MAEが最小となるのはが
の中央値になったときです。そのため、MAEを最小にしようとすることは予測値を中央値に近づけることに繋がります。
「MAEが最小となるのはが
の中央値になったとき」の証明については、以下で検討しております。（ちゃんとした証明はわからなかったので雰囲気を記載しております。）
RMSEは予測誤差の2乗の平均の２乗根です。ルートを取らなければMean Squared Error (MSE)となります。これを評価指標としてモデル選択を行うと、その予測モデルは平均値を予測するモデルとなります。これは「RMSEを評価指標として使う」= 「RMSEが最小となるモデルを選択する」ということになるのですが、RMSEが最小になるのはが
の平均値になったときです。そのため、RMSEを最小にしようとすることは予測値を平均値に近づけることに繋がります。（これはMSEを評価指標にした時も同じです。）
「RMSE(or MSE)が最小となるのはが
の平均値になったとき」の証明については、以下で検討しております。
平均値がMSE(or RMSE)を最小化する理由 - はだだだだ
評価指標自体はMAEの方がわかりやすいと思いますが、導かれる予測値の性質としてはRMSEの方がわかりやすいと思います。MAEもRMSEも元の予測対象となる時系列の単位をそのまま引き継いでいるため、異なる時系列間でMAEやRMSEを比較することはできないです。ただし、異なる時系列間でMAEやRMSEを比較するケースは、予測手法の性能や特徴を評価する時が主なケースかな？と思いますので、実際に予測モデルを作るとき（ある１つの時系列データの予測モデルを作るとき）は関係ないかと思います。
個人的には実務等で予測モデルを作るときはとりあえずMAEかRMSEを使用すれば良いのではないかと思います。
Percentage Errors
Percentage Errorsの評価指標の例とし、以下が挙げられます。
- Mean Absolute Percentage Error : MAPE =
- symmetric Mean Absolute Percentage Error : sMAPE =
MAPEは予測誤差を実測値で割って%表記にしたものです。これを使用することで異なる時系列データに適用した予測モデルの性能を比較することができます。ただし、分母となる実測値が0の場合は割り算が定義できず、また、実測値が0に近い小さな値の場合、MAPEはとても大きな値になるため、評価指標として不安定になります。
また、MAPEを使用する際の注意点として、そもそも予測対象となる時系列データが「%表記して問題ないか？」ということが挙げられます。例えば、気温のデータについて、「東京の天気が昨日より10%上がりました」などと言うことは無いと思います。これは、気温というもの区切り（1℃と2℃）が恣意的で、必ずしも2℃が1℃の2倍暑いとは言えないためです。実際、華氏で表すと2℃は35.6度、1℃は33.8度のため、2倍では無くなります。
そのため、温度を予測するモデルを作る時には、評価指標としてMAPEは使用できません。
MAPEについてはもう1点注意すべき点があります。例えば、,
のとき、
になるのに対し、
,
のとき、
になります。このように実測値と予測値が入れ替わっただけで（上回ったか、下回ったかの違いだけで）、評価指標の値が変わってしまいます。
この点をケアするために考案されたのが、 SMAPEです。先ほどの例では、実測値と予測値の大小関係は関係なく、いずれもになります。
ただし、sMAPEには別の問題があり、以下の記事と論文では推奨されていません。
Another look at measures of forecast accuracy | Rob J Hyndman What the MAPE is FALSELY blamed for, its TRUE weaknesses and BETTER alternatives! | R-bloggers
一方、sMAPEはKDD CUP 2018で評価指標として使われていたりするため、ポピュラーな指標ではあるようです。
KDD CUP of Fresh Air - Biendata
Scaled Errors
最後にScaled Errorsです。これはPercentage Errorsと同様、異なる時系列データの予測モデルの性能を比べる際に使用できます。Scaled Errosの例として以下が挙げられます。
- Mean Absolute Scaled Error : MASE =
ただし
MASEは予測モデルの予測誤差と、ベンチマークとなるナイーブ予測の予測誤差を比較し、その比率を取っています。季節性が無い時系列の場合のナイーブ予測として、t-1の値をそのままtの予測値として使用する場合が上記のの”季節性が無い場合”に当たります。（ナイーブ予測の方法は任意？だと思います。）
こうすることで、MASEが1より小さければ予測モデルはナイーブ予測に勝っている（予測誤差が小さい）とわかります。この指標では分母にナイーブ予測の予測誤差を使用することで単位の影響が無くなり、異なる時系列データの予測モデルの性能比較に使用できます。
なお、MASEを最小化するモデルはMAEと同様、中央値を予測するようなモデルになります。
Another look at measures of forecast accuracy - ScienceDirect
自分なりの結論
実務で予測モデルを作る際には、シンプルでわかりやすいScale-dependent Errorsを使用すれば良いと思います。MAEとRMSEのどちらを使うかについては、「両方試してみて、両方良さそうなモデルを選択する」というやり方が、現実的な方法ではないかと思います。
また、複数時系列データへの性能評価を行いたい場合はMASEを使用すれば良いのではないかと思います。
補足
時系列モデルの予測精度を競う大会ではモデル評価指標として何を使うかが重要になります。こちらの記事では時系列モデルの予測精度を競う大会の歴史について記述されており、モデル評価指標の歴史についてもわかりやすくまとめられています。
平均値がMSE(or RMSE)を最小化する理由
厳密な証明かはわかりませんが、以下に自分なりの理解をまとめます。MSEの最小化問題はRMSEの最小化問題と解が同じ(?)だと思いますので、MSEの方で検討します。
(記法)
:
のベクトル
:
のベクトル
:
のベクトル
ここで(5)の第3項に注目して
となるので、
となります。
ここで(6)をに関する式とみなすと、
が含まれるのは第2項のみとなるため、第2項を最小化する
がMSE(or RMSE)を最小化する
となります。
第2項は2乗された値のため、最小値は0であり、第2項が0になるのは、]のときです。
]は
の期待値（平均値）のため、「平均値がMSE(or RMSE)を最小化する」と言えます。
(参考資料)
* https://scholar.harvard.edu/files/danielyewmaolim/files/api-208section1.pdf
中央値がMAEを最小化する理由
厳密な証明はよくわからなかったので、自分なりのざっくりした理解を記載します。
(1)(2)予測誤差を式変形していきます。
(3)が連続で
の値をとる場合、絶対値の期待値は積分で表せます。このとき、
を確率変数と考え、
は定数とみなすことにします。
は
の確率密度関数です。
(4)積分区間をで分けます。
(5)大小関係に気をつけて絶対値を外します。
(6)MAEを最小化する予測値を考えます。これまで
は定数とみなしていましたが、ここからは変数とします。MAEを
で最小化する場合、1階の微分条件はこの通りになります。
(7)微分を行います。
(8)移項します。
(9)(8)の式を確率の記号で読み変えます。
(9)のように表せるとき、は
の中央値となります。これは連続変数の場合の中央値の一般的な定義です。
以上からが
の中央値の時、MAEが最小化されると言えそうです。（中央値がMAEを最小化すると言えそうです。）
参考資料
ホームページをレンタルサーバーからGitHub Pagesに移してみた
これまでXFREEのレンタルサーバーにホームページを作り、そこにRの勉強で作った資料などをあげていました。無料で広告も入らないため不満はなかったのですが、3ヶ月に1度更新処理（ログインして更新ボタンを押す）が必要なことと、更新するたびにローカルでhtmlファイル作成→FTPソフトでファイルのアップロードを行う必要があるため、少し面倒だなと思っておりました。
そこで、GitHub Pagesにホームページを移動することにしました。
また、ホームページの内容をＲで管理できるようにR Markdownでホームページを書くことにしました。
移動前のホームページ
http://hadadada00.html.xdomain.jp/
（2019/11/30以降は期限切れで見られなくなると思います。）
見た目はこんな感じです。
XFREEサーバーの中身はこんな感じです。
index.htmlはこんな感じです。
<!DOCTYPE html> <html lang="ja"> <head> <meta charset="utf-8"> <title>hadadada00</title> </head> <body> <h1>hadadada00</h1> <h3>blog</h3> <a href="https://hadadada00.hatenablog.com/">はだだだだ</a><br> <h3>report</h3> <a href="http://hadadada00.html.xdomain.jp/gsw_2016_2017.nb.html">Data analaysis on Golden State Worriors in 2016-2017 season.</a><br> <a href="http://hadadada00.html.xdomain.jp/NBA_players_height.html">NBA playears' height.</a><br> <a href="http://hadadada00.html.xdomain.jp/homestates.html">NBA playears' home states.</a><br> <h3>data</h3> <a href="http://hadadada00.html.xdomain.jp/gdp.csv">gdp.csv</a><br> <h3>web application</h3> <a href="https://hadadada00.shinyapps.io/google_trends/">Google Trends downloader</a><br> <br> <br> <br> <footer>(c) 2019 hadadada00</footer> </body> </html>
レポジトリの作成
githubにログインしてレポジトリの新規作成を行います。
ディレクトリ名を[username].github.ioにします。
ホームページ用のR Projectの作成
次に、GitHub Pagesの内容をＲで管理できるように、R Projectを作成します。
R Studio > New Project > Version Cotrol > Git と進み、先ほど作成したリポジトリのURLを入力します。
これでR studioからホームページを更新できるようになりました。
ホームページをR Markdownで作成する。
R Markdownファイルを格納するためのRmdフォルダを作成します。
Rmdフォルダ内に、index.Rmdという名前でファイルを作り、このスクリプトからホームページのトップページになるindex.htmlを作成します。
index.Rmd
--- title: "hadadada00" author: "hadadada00" date: "2019年9月1日" output: html_document knit: (function(inputFile, encoding) { rmarkdown::render(inputFile, encoding = encoding, output_dir = "../") }) --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` ## blog - [はだだだだ](https://hadadada00.hatenablog.com) ## report - [Data analaysis on Golden State Worriors in 2016-2017 season](https://hadadada00.github.io/html/gsw_2016_2017.nb.html) - [NBA playears' height](https://hadadada00.github.io/html/NBA_players_height.html) - [NBA playears' home states](https://hadadada00.github.io/html/homestates.html) ## data - [gdp.csv](https://hadadada00.github.io/csv/gdp.csv) ## package - [nbastats](https://github.com/hadadada00/nbastats) ## web application - [Google Trends Downloader](https://hadadada00.shinyapps.io/google_trends/)
ポイントはyamlヘッダに以下を追加していることです。
knit: (function(inputFile, encoding) { rmarkdown::render(inputFile, encoding = encoding, output_dir = "../") })
参考：
r - Rmarkdown directing output file into a directory - Stack Overflow
これをすることでR Markdownファイルは/hadadada00.github.io/Rmdフォルダの中に格納しておき、生成されるindex.htmlは1つ上の/hadadada00.github.ioフォルダに作成されます。
後は以前XFREEサーバーにアップしていたhtmlファイルやcsvファイルをそれぞれhtmlフォルダとcsvフォルダの中にコピーしてRStudioからgithubにcommit > pushをすれば完了です。
ローカルのR Projectのフォルダ構成は以下です。
実際に作成したホームページが以下です。
https://hadadada00.github.io/
以上
RでARモデルの勉強
時系列分析の勉強をはじめることにしました。まずはARモデルです。
以下の記事を見ながらARモデルの構築と、モデル検証の勉強をしました。
データ
Google Trendsのヒット数を使用します。検索ワードは”温泉”にしました。
期間は2016年1月1日～2019年1月1日のdailyデータで、データの取得には、以下で作成した自作ツールを使用しました。
hadadada00.hatenablog.com
基礎分析
取得したcsvファイルをRに読み込んで簡単にトレンドと季節性がないか確認しました。
最初にトレンドを確認します。
# import pacages ---------------------------------------------------------- library(lubridate) library(tidyverse) # read -------------------------------------------------------------------- onsen <- read_csv("./onsen_20160101-20190101.csv") %>% mutate(date = as.Date(date)) # check data -------------------------------------------------------------- # 1. trend onsen %>% ggplot(aes(x = date, y = hits_adj)) + geom_line() + scale_x_date(date_labels = "%Y%m%d")
全体的に少し右上がりか横ばいです。
次に曜日による季節性を確認しました。
# 2. weekly seasonality onsen %>% mutate(week_day = wday(date, label = TRUE)) %>% ggplot(aes(x = week_day, y = hits_adj)) + stat_summary(fun.y = "mean", color = "red", geom = "point") + stat_summary(fun.y = "median", color = "blue", geom = "point")
明らかに曜日による違いがあります。日曜日と土曜日が多く、次に月曜と金曜が少し多めで、火曜、水曜、木曜が少ないです。
温泉というレジャーに関する検索ヒット数のため、休日に数が増えるのは妥当だと思います。
最後に月による季節性を確認します。
# 3. monthly seasonality onsen %>% mutate(month = month(date, label = TRUE)) %>% ggplot(aes(x = month, y = hits_adj)) + stat_summary(fun.y = "mean", color = "red", geom = "point") + stat_summary(fun.y = "median", color = "blue", geom = "point")
こちらも月による季節性がありそうです。１月、８月が多いのは冬休みと夏休みがあり温泉旅行に関心が高まるからだと思います。（夏でも多いのが少し意外です。）それ以外は概ね暑い時期は少なく、寒い時期は多いという傾向が見られます。
モデル構築
まずはAR(p)のp(次数)の目安を得るために、自己相関係数と偏自己相関係数を見ます。
stats::acfとstats::pacf関数をそれぞれ使用します。
自己相関係数を見ると、自己相関が存在し、特に７日前、１４日前、と１週間区切りで相関が高いことが確認されます。
次に偏自己相関係数を確認します。AR(p)モデルの次数はこちらの結果が参考になります。
（例えば、自己相関係数だとt=1とt=7の相関を見るときにt=2, t=3など他の時点の影響が含まれてしまいます。本当はt=1, t=7間に直接相関がなかったとしてもt=1, 2, 3, 4, 5, 6, 7と影響が連鎖していき、見た目上t=1とt=7に相関があるように見えてしまいます。そのため次数の参考には偏自己相関係数を使用します。）
結果を見るとt=-1に大きなプラスの相関があります。
以降はt=-6, -12, -13, -20と概ね6～8日おきに正の相関がある一方、t=-8, -15, -22は負の相関がでております。
前日と、概ね1週間単位のラグに影響を受けていることがわかります。
モデリング
stats::ar関数でモデルを構築します。デフォルトでは次数がAICを基準に自動的に選ばれます。引数としては最大次数のみ個別に設定します。
デフォルトでは以下の記述にあるように、min(N-1, 10*log10(N) N=number of observations の値がセットされるため、2016年1月1日～2019年1月1日のdailyデータの場合、N=1097で、デフォルトの最大次数は30になります。
order.max
maximum order (or order) of model to fit. Defaults to the smaller of N-1 and 10*log10(N) where N is the number of non-missing observations except for method = "mle" where it is the minimum of this quantity and 12.
基礎分析で月の季節性があることが確認できたため、最大1年分の影響を検討対象にするべく、最大次数は400としました。
# model building ---------------------------------------------------------- # AR model ar_model <- ar(onsen$hits_adj, order.max = 400) summary(ar_model) > summary(ar_model) Length Class Mode order 1 -none- numeric ar 30 -none- numeric var.pred 1 -none- numeric x.mean 1 -none- numeric aic 401 -none- numeric n.used 1 -none- numeric n.obs 1 -none- numeric order.max 1 -none- numeric partialacf 400 -none- numeric resid 1097 -none- numeric method 1 -none- character series 1 -none- character frequency 1 -none- numeric call 3 -none- call asy.var.coef 900 -none- numeric
summary()の実行結果のうち、arが次数を表わしています。arの値は30のため、AR(30)のモデルが選択されました。
モデル検証
以下の2点を確認します。
- 誤差項の独立性
- 誤差項の正規性
これらはいずれもAR(p)モデルの前提になっている条件です。
まずは誤差項の独立性です。Ljung-Box検定を行います。
# check model ------------------------------------------------------------- # 1. Ljung-Box test(1978) (independency of error term) # H0 : data doesn't show auto correlation # small p-value means there is a autocorrelation Box.test(ar_model$res, type = "Ljung") >Box.test(ar_model$res, type = "Ljung") Box-Ljung test data: ar_model$res X-squared = 3.7343, df = 1, p-value = 0.05331
Ljung-Box検定の帰無仮説は「誤差項に自己相関が存在しない」です。
p値はは0.053のため、5%有意水準で帰無仮説は棄却されませんが、10%有意水準では棄却されます。そのため、検定結果としては誤差項に自己相関がないとはいえないです。
次に誤差項が正規分布に従うか検証します。Jauque-Bera検定を行います。
このとき、入力データに欠損がでないよう気をつけます。
# 2. Jauque-Bera test(1987) (normality of error term) # H0 : data has a normal distribution # small p-value means the data doen't follow normal distribution # before testing model, remove NA data in dataset jarque.bera.test(ar_model$res %>% tail(n = -30)) >jarque.bera.test(ar_model$res %>% tail(n = -30)) Jarque Bera Test data: ar_model$res %>% tail(n = -30) X-squared = 1533.9, df = 2, p-value < 2.2e-16
p値は小さいため、有意水準1%以下で棄却され、誤差項は正規分布に従わないという結果がでました。
検証結果を見ると、誤差項にまだ自己相関がのこっており、正規分布にも従っていないため、モデリングにはまだ足りない要素があるといえます。
今回はモデリングの手順確認が主目的のため、追加のモデル構築は行いません。
予測
構築された（イマイチな）モデルを使用して、予測を行ってみます。
ただし今回はモデルの精度を検証するというよりはARモデルで予測するとどのような予測値が得られ、どのような限界があるかを確かめることに重点をおきます。
そのため以下のようなやり方でできるだけ長期間の予測値を得られるようにします。
- 2016-01-01～2016-01-30の実際の値(hits_adj)を使用して以降の期間の予測（hits_hat）を作成する。
- 2016-01-31の値は2016-01-01～2016-01-30の実測値を使って予測する。
- 2016-02-01の値は2016-01-02～2016-01-30の実測値と2016-01-31の予測値を使用して予測する。
- 以降同様に予測に使用するデータを実測値から予測値にずらしていく。
# predict ----------------------------------------------------------------- input <- head(onsen$hits_adj, n = 30) pred <- predict(ar_model, newdata = input, n.ahead = (length(onsen$hits_adj) - 30)) onsen_hits_hat <- union_all(rep(NA, 30), pred$pred) # data frame for plot df <- data.frame(date = onsen$date, hits_adj = onsen$hits_adj, hits_hat = onsen_hits_hat) # plot true value and predicted value df %>% ggplot() + geom_line(aes(x = date, y = hits_adj)) + geom_line(aes(x = date, y = hits_hat), color = "red") + xlim(as.Date("2016-01-01"), as.Date("2019-01-01"))
実際の値と予測値をプロットすると以下のようになります。
と予測値の比較（３年分）
まず3年後は明らかに水準も違いますし、上下の変動幅も全く異なります。
描画期間を変更して、短期間での予測値の動きを見てみます。
# plot true value and predicted value df %>% ggplot() + geom_line(aes(x = date, y = hits_adj)) + geom_line(aes(x = date, y = hits_hat), color = "red") + xlim(as.Date("2016-01-01"), as.Date("2019-07-01"))
描画期間を2016-01-01～2016-07-01に変更しています。予測値は2016-01-31からのため、概ね半年間の予測をしています。
結果をみると、最初の3ヶ月あたりは、上昇の動きと下落の動きを捉えられており、かつピークとボトムの水準の誤差はそれほど大きくないため、あれば予測値としてある程度は使えると思います。
しかし、3ヶ月を過ぎると、上昇、下落のタイミングは捉えられているものの、全体的な水準が異なっているため、予測値としては使えません。
今回構築されたモデルがAR(30)のため、曜日の季節性はとらえられているものの、月の季節性が捉えられていないことが原因と思われます。
とはいえ、必要な予測期間が3ヶ月程度であればARモデルというシンプルなモデルでも、それなりに使える予測値が得られるということがわかりました。
shinyで簡単なwebアプリを作ってshinyapp.ioに公開してみた
shinyの練習をしました。以下の記事で作成したコードを流用して、google trendsのdailyの結果をcsv形式で取得するウェブアプリを作ります。
具体的には以下のことができるようにしました。
- 検索ワードを指定
- 期間を指定
- 実行ボタンを押すと、条件に沿ったデータを取得
- 取得したデータをグラフで表示（確認用）
- ダウンロードボタンを押すとcsv形式でデータがダウンロードされる
参考記事
以下の記事が大いに参考になりました。(公式チュートリアルの日本語訳)
RStudio Shiny チュートリアル レッスン1 ようこそ Shiny へ - Qiita
app.Rファイルの作成
# import packages --------------------------------------------------------- library(gtrendsR) library(shiny) library(tidyverse) # ui ---------------------------------------------------------------------- ui <- fluidPage( # title titlePanel("Google Trends downloader"), # Sidebar with a slider input for number of bins sidebarLayout( sidebarPanel( # keyword textInput("keyword", label = h5("Keyword"), value = ""), # start date textInput("start", label = h5("Start date (yyyy-mm-01)"), value = ""), # end date textInput("end", label = h5("End date (yyyy-mm-01)"), value = ""), # action button actionButton("update", "Get data") ), # Show a plot of the generated distribution mainPanel( # plot line plotOutput("linePlot"), # Download Button downloadButton("downloadData", "Download") ) ) ) # server ------------------------------------------------------------------ server <- function(input, output) { df <- eventReactive(input$update, { # prepare arguments for grends() function start_date <- seq(as.Date(input$start), as.Date(input$end), by = "month") end_date <- seq(as.Date(input$start), length = length(start_date), by = "month") period <- cbind(data.frame(s = start_date[-length(start_date)]), data.frame(e = end_date[-1])) %>% mutate(period = str_c(s, " ", e)) %>% select(period) # get google trends data trends <- lapply(period$period, gtrends, keyword = input$keyword, geo = "JP") # extract index data trends_df <- trends %>% map("interest_over_time") %>% bind_rows() %>% mutate(weight_temp = coalesce(ifelse(date == lag(date, 1), lag(hits, 1) / hits , 1), 1), weight = cumprod(weight_temp), hits_adj = hits * weight) %>% distinct(date, .keep_all = TRUE) %>% select(date, hits_adj) # output trends_df }, ignoreNULL = FALSE) output$linePlot <- renderPlot({ req(input$update) df() %>% mutate(date = as.Date(date)) %>% ggplot(aes(x = date, y = hits_adj)) + geom_line() + scale_x_date(date_labels = "%Y%m%d") }) output$downloadData <- downloadHandler( filename = function() { paste(input$keyword, ".csv", sep = "") }, content = function(file) { write_csv(df(), file) } ) } # run --------------------------------------------------------------------- shinyApp(ui = ui, server = server)
アプリの実行
まず、Rstudio上で"Run app"を押してアプリを実行します。
すると、別ウィンドウにアプリが表示されます。
検索ワードと期間をしていして"Get data"を押します。
するとgoogle trendsからデータが取得され、グラフが描画されます。
最後にDownloadボタンを押すと、csvファイルがダウンロードされます。
デプロイ
今回はshinyapp.ioにデプロイすることにしました。
やり方はいくつかあるようですが、Rstudioから簡単にできました。
アプリを実行すると右上にpublishボタンがでてくるので押します。（直前の画像を参照）
押すとshinyapp.ioと繋げるためのダイアログボックスがでてきます。
（なお、事前にshinyapps.ioに登録してtokensを取得しておいてください。）
Nextを押します。
ShinyApps.ioを選択します。
shinyapps.ioのアカウントページで発行されるtokensのスニペットを貼り付けます。
アカウントが紐付きました。publishを押せばあとは自動的にデプロイしてくれます。
デプロイ後はwebで公開されるため、以下のように自由にアクセスできます。
Google Trends downloader
以上
google trendsで長期間のdailyデータを取得する方法
google trendsのデータを使って時系列分析の勉強をしようと思いました。gtrendsRというRパッケージを使用すればRからAPI経由でデータを取得できるようです。
gtrendsRの関連記事
https://mrunadon.github.io/gtrendsR/
しかし、google trendsの仕様ではdailyで取得できるのは最大9ヶ月程度で、かつ、数値は抽出時の最大値を100として指数化されているため、期間を分けて抽出したとしても単純に比較が出来ません。
そこで、以下の内容を実現するための関数を別途作成しました。
- startとendで期間を指定可能（9ヶ月以上可）
- 抽出した数値の基準が揃っており、長期間で比較が可能
コードの全体は以下です。
# import packages --------------------------------------------------------- library(gtrendsR) library(lubridate) library(tidyverse) # parameters -------------------------------------------------------------- start <- "2018-04-01" # yyyymm (only 1st date of month) end <- "2019-08-01" # yyyymm (only 1st date of month) keyword <- "HKT48" geo <- "JP" # get google trends data -------------------------------------------------- # change parameters into aruguments of gtrends() function start_date <- seq(as.Date(start), as.Date(end), by = "month") end_date <- seq(as.Date(start), length = length(start_date), by = "month") period <- cbind(data.frame(s = start_date[-length(start_date)]), data.frame(e = end_date[-1])) %>% mutate(period = str_c(s, " ", e)) %>% select(period) # get google trends data trends <- lapply(period$period, gtrends, keyword = "usj", geo = "JP") # extract index data trends_df <- trends %>% map("interest_over_time") %>% bind_rows() %>% mutate(weight_temp = coalesce(ifelse(date == lag(date, 1), lag(hits, 1) / hits , 1), 1), weight = cumprod(weight_temp), hits_adj = hits * weight) %>% distinct(date, .keep_all = TRUE) trends_df %>% mutate(date = as.Date(date)) %>% ggplot() + geom_line(aes(x = date, y = hits_adj)) + scale_x_date(date_labels = "%Y%m%d")
処理の流れ
期間してについてはyyyy-mm-01形式しか対応していませんが、機関を余分に設定すれば欲しい期間のデータはそろうため、このままにしております。
流れは以下です。
- パラメータで検索ワードと期間(start, end)を指定
- startとendから1ヶ月区切りのパラメータベクトルperiodを作成
- periodの分だけgtrendsR::gtrendsをまわす
- 取得した複数のデータフレームを連結し、検索数(hits)を調整したhits_adjを作成
処理のポイント
1ヶ月区切りの日付を作成する歳に、1日だけ日付が重なるようにしています。
# parameters -------------------------------------------------------------- start <- "2018-04-01" # yyyymm (only 1st date of month) end <- "2019-08-01" # yyyymm (only 1st date of month) keyword <- "HKT48" geo <- "JP" # get google trends data -------------------------------------------------- # change parameters into aruguments of gtrends() function start_date <- seq(as.Date(start), as.Date(end), by = "month") end_date <- seq(as.Date(start), length = length(start_date), by = "month") period <- cbind(data.frame(s = start_date[-length(start_date)]), data.frame(e = end_date[-1])) %>% mutate(period = str_c(s, " ", e)) %>% select(period) >period > period period 1 2018-04-01 2018-05-01 2 2018-05-01 2018-06-01 3 2018-06-01 2018-07-01 4 2018-07-01 2018-08-01 5 2018-08-01 2018-09-01 6 2018-09-01 2018-10-01 7 2018-10-01 2018-11-01 8 2018-11-01 2018-12-01 9 2018-12-01 2019-01-01 10 2019-01-01 2019-02-01 11 2019-02-01 2019-03-01 12 2019-03-01 2019-04-01 13 2019-04-01 2019-05-01 14 2019-05-01 2019-06-01 15 2019-06-01 2019-07-01 16 2019-07-01 2019-08-01
これによって、各期間で取得したデータフレームの検索数（hits）の調整ができるようになります。
調整をおこなっているのは以下の箇所です。
# get google trends data trends <- lapply(period$period, gtrends, keyword = "usj", geo = "JP") # extract index data trends_df <- trends %>% map("interest_over_time") %>% bind_rows() %>% <b> mutate(weight_temp = coalesce(ifelse(date == lag(date, 1), lag(hits, 1) / hits , 1), 1), weight = cumprod(weight_temp), hits_adj = hits * weight) </b>
まず、取得したデータフレームを縦につなげます。
trends_df <- trends %>% map("interest_over_time") %>% bind_rows() %>%
次に、日付が連続しているところ（元々データフレームの境目）のhitsの水準を私用して、前後のデータフレームの調整係数を作成します。
mutate(weight_temp = coalesce(ifelse(date == lag(date, 1), lag(hits, 1) / hits , 1), 1), weight = cumprod(weight_temp),
ここの処理を表にすると以下のようになります。
日付が連続しているところでweight_tempを作成し、その連乗積を計算することでweightを計算しています。
そしてhitsにweightをかけることで最初に取得したデータフレーム内のhitsの水準に、以降のデータフレームのhitsを合わせています。
実行例
以下のように、長期のdailyのデータが取得できています。
以上
