R言語(プログラミング)を用いて非線形成長モデルの将来予測・考察をしてみた(エストニアの政府データを利用）
Tere!🇪🇪
北欧はエストニアの、タルトゥ大学で計量経済学を学んでいるとみたです。
４月になり、エストニアは急に暖かくなってきました。
つい最近までコートで出かけていたのに、今じゃ半袖の人もちらほら見かけます。
ここ留学が決まった時は、色んなことをブログに書こうと思ったのですが、全然書けていませんね....（シンプル怠惰）
来月あたり、"エストニア留学生の１日”と題して、それっぽい動画をあげようかな...とは思っています。
（先に宣言して自分を追い込むタイプ。故に有言不実行が多い。良くない。）
ということで今回も、主に自分の備忘録として、ニッチなプログラミングと経済学に関する記事を書いていきます〜。
前回に引き続き、エクセルでも出来る事なのですが、
モデル作成から描画までR言語を用いるととてもスピーディーなので、個人的に良いなと思っています。
どっかのRを使われされている大学院生とかに役立てばいいなあ。
前置き
今回のテーマは、エストニアの政府データを利用して、X-Roadの経済成長のForecasting（予測）を行なってみよう！といったところです。
X-Roadというのは、エストニアの電子政府を根源で支えているとっっても重要なサービスです。
eID(日本でいうマイナンバー）にくくりつけられた、運転免許や処方箋と言った個人の様々な情報が、X-Roadという基盤を通して、必要とする機関に開示されます。その結果、個人は運転免許証やお薬手帳を持ち歩く必要がありませんし、納税もパソコンでサクッと行うことができます。
そんなことしたら、個人情報の安全性はどうなるの〜〜とお思いになるかもしれませんが、ここで肝なのが、X-Roadは、個人の情報を一括して保持しているのではなく、各機関がその情報を必要とし、情報開示の要求を発行した際に、情報を得る権利を各機関に与えるという役割を持っているという点です。これによって、情報は安全に透明性の高いシステムの中で運用されます。
（もっとうまく説明されている記事がたくさんあるので、気になった方はググってみてくださいね。）
説明が長引いて、ちょっと脱線してしまいましたね（汗）
今回は、そんなX-Roadの経済成長に関して、R言語を用いて予測・考察をしていきます。
目的
具体的に何を行うかというと、X-Roadに関する２種類の予測と考察をしていきます。
経済予想１（情報開示請求数）
データ公開時までの非線形成長モデルが維持された場合、例えば2020年12月にX-Roadにおいて行われる情報開示請求数はいくつになるだろうか。(データから成長モデルも求める)
経済予想２（サービス数）
データ公開時までの線形成長モデルが維持された場合、例えば2019年3月にX-Roadのサービス（例：医師と薬剤師間の患者の処方箋の共有サービス）の数のはいくつになるだろうか。(データから成長モデルも求める)
予測のプロセス
どちらの予測のプロセスも、大まかに以下の３つフェーズに分類できます。それぞれRを用いてスピーディーに行います。
①標本数と時間をインプットしプロットを作成
②成長モデルを求める（パラメーターを推定）
③予測を行う
2020年12月次の情報開示請求数の予測
標本数と時間をインプットしプロットを作成
まずは、データをインプットします。関数option()のscipenによって指数表記を回避しておきましょう。
options(scipen=100000) xroad=read.table("http://www.ut.ee/kristjan.vassil/wp-content/Xroad_monthly.csv",,header = T,sep=",")
そうしたら、プロットする前に、表（ベクトル）を定義しましょう。
y軸を請求数、x軸をデータが始まっている2003年からの月数とします。
# Define vectors y (queries) and x (count of months from Jan, 2003) y=xroad$Queries x=1:length(y)
定義したら、プロットします。
plot(x,y,ylab="Queries",xlab="Month")
するとこんな風に描かれます。
成長モデルを求める（パラメーターを推定）
次に、の成長モデルを求めていきます。
#Estimate parameters a and b for y=ab^x exponent.queries=nls(y~a*b**x,start=list(b=1,a=21670)) exponent.queries
するとこんな風に係数が求まります。
Nonlinear regression model model: y ~ a * b^x data: parent.frame() b a 1.023 1659219.979 residual sum-of-squares: 1929083179148436 Number of iterations to convergence: 11 Achieved convergence tolerance: 0.0000009646
そうしたら、このモデルを描画させましょう。
#Plot fitted line plot(x,y,ylab="Queries",xlab="Month") lines(x,coef(exponent.queries)["a"]*coef(exponent.queries)["b"]**x,lwd=2)
すると、こんな図ができます
予測を行う
では、最後に予測を行なっていきます。
months2=1:216 # 予測をする分の216ヶ月目まで範囲を延ばす plot(months2,coef(exponent.queries)["a"]*coef(exponent.queries)["b"]**months2,type="l",ylab="Queries",xlab="Month",lwd=2,lty=2) # draw a line using the above sepcified function with estimated parameters points(x,xroad$Queries) # add actual data lines(x,coef(exponent.queries)["a"]*coef(exponent.queries)["b"]**x,lwd=2) # mark prediction line with dashed line queries216=(coef(exponent.queries)["a"]*coef(exponent.queries)["b"]**months2) # calculate predicted values for all months sum(queries216[months2>=216]) # 216ヶ月目の月間合計予測データを計算 ||< [f:id:tommy-acoustic-7:20190427075034j:plain] こんな予測グラフが出力されました。 結果としては、 >|r| [1] 215047461
より、データ公開時までの非線形成長モデルが維持された場合、2020年12月にX-Roadにおいて行われる情報開示請求数は215047461になることがわかった。
2019年3月次のサービス数の予測
ほとんど上記と同じ流れで行なっていきます。違うのは線形か非線形かの違いくらい。
標本数と時間をインプットしプロットを作成
すでにインプットは行われているので、いきなりプロットを描画しちゃいます。
#Step 1 - plot the relationship plot(x,xroad$Services,ylab="Services",xlab="Month") #Plot services over time
こんな感じですね。
成長モデルを求める（パラメーターを推定）
そしたら、線形モデルを求めましょう。
#Step 2 - Estimate parameters a and b for y=a+bx linear.services=lm(xroad$Services~x) linear.services
すると係数が求まります。
Call: lm(formula = xroad$Services ~ x) Coefficients: (Intercept) x 55.80 10.27 ||< では、グラフを出力していきます！ >|r| #Plot fitted line plot(x,xroad$Services,ylab="Services",xlab="Month") abline(linear.services,lwd=2)
予測を行う
実際、手作業でもできちゃうレベルなのですが、まあここは基礎を学ぶ意味も込めてR言語で予測していきます。
months2=1:195 plot(months2,predict(linear.services,newdata=data.frame(x=months2)),type="l",ylab="Services",xlab="Month",lwd=2,lty=2) # draw a line using the above sepcified function with estimated parameters points(x,xroad$Services) # add actual data lines(x,predict(linear.services,newdata=data.frame(x=x)),lwd=2) # mark prediction line with dashed line linear.services=55.79820+10.26877*months2 # calculate predicted values for all months linear.services[181:195] ||< このコードによって、181ヶ月目から195ヶ月目までのy切片（サービス数）の値を出力してくれますね。 >|r| [1] 1914.446 1924.714 1934.983 1945.252 1955.521 1965.789 1976.058 1986.327 1996.596 2006.864 2017.133 2027.402 2037.671 2047.940 2058.208
したがって、2058.208個のサービスが利用されるという予測が立ちました。
図はこんな感じになります。
考察（消費者数）
ここまではX-Roadの成長の予測を行なってきたのですが、ここではついでに消費者数のプロットを打って、得られる図から考察を行います。
既にインプットしてあるデータから、消費者数の推移を描画させます。
このコードから、
plot(x,xroad$Consumers,ylab="Consumers",xlab="Month")
明らかに線形モデルだとフィットしませんね。
消費者の増加がこのような曲線になる理由は、イノベーター理論が説明可能でしょう。
イノベーター理論とは、1962年に米・スタンフォード大学の社会学者、エベレット・Ｍ・ロジャース教授（Everett M. Rogers）が提唱したイノベーション普及に関する理論で、新しいサービスがどのように消費者の間で普及していくかという説明をするものです。
以下のグラフのように、最初は新しいもの好きのイノベーティブな消費者からゆっくり広まり、みんながやっているならやるか、という人たちの間でグッと広がり、最後にまあシャーなしでやるかという風に変化をあまり好かない人たちの間で広まっていきます。
現実に、社会科学の理論が成り立っているのを観察することのできる面白い一例ですね。
