www.crosshyou.info

政府統計の総合窓口のデータや、OECDやUCIやのデータを使って、Rの練習をしています。ときどき、読書記録も載せています。

都道府県別の老人福祉費と児童福祉費の分析1 - R言語で都道府県別の平均値を算出する。

今回は、都道府県別の老人福祉費と児童福祉費の分析をしてみようと思います。

データは、政府統計の総合窓口e-stat(www.e-stat.go.jp)から取得しました。

f:id:cross_hyou:20190926190817j:plain

取得したデータは5つです。総人口、総面積、県内総生産額、老人福祉費、児童福祉費です。

f:id:cross_hyou:20190926190908j:plain

ExcelのCSVファイルにこんな感じでダウンロードしました。

f:id:cross_hyou:20190926190946j:plain

R言語のread.csv関数でデータファイルを読込み、str関数で構造を確認します。無事に読み込んでいます。

まず、NAの行をomit.na関数で削除します。

f:id:cross_hyou:20190926191250j:plain

summary関数で平均値などを表示しています。NAが無くなっています。

年度と地域で全部そろっているか確認しましょう。table関数を使います。

f:id:cross_hyou:20190926191646j:plain

table関数を使う前にas.characcter関数とas.vector関数で年度の不要なファクタを削除しています。2006年度から2015年度の10年間ですべて47となっていますので47都道府県全部のデータがそろっています。

table関数をPrefでもやってみましょう。

f:id:cross_hyou:20190926192017j:plain

47都道府県すべて10個のデータがありますから、データはそろっていることが確認できました。

過去10年間平均の老人福祉費を都道府県別に調べましょう。tapply関数を使います。

f:id:cross_hyou:20190926192526j:plain

tapplyで都道府県別の平均値を計算し、avgOldというオブジェクトに格納しています。そして、sort関数で値の小さい順に表示しています。鳥取県が一番すくなく、169億8846万1千円です。東京都が2635億5029万6千円です。

児童福祉も同じようにtapply関数で調べましょう。

f:id:cross_hyou:20190926193419j:plain

香川県が児童福祉費は一番少ないですね。97億9266万7千円です。

東京都が一番多く、1980億1086万9千円です。児童福祉費よりも老人福祉費が多いのですね。

ついでに、総人口、総面積、県内総生産の10年間の平均値も出しておきましょう。

f:id:cross_hyou:20190926193929j:plain

総人口の一番少ないのは鳥取県です。58万8110人です。一番多いのは東京都で1313万8565人です。

f:id:cross_hyou:20190926194235j:plain

香川県が一番小さいです。18万7656haです。一番大きいのは北海道で、784万2076haです。

f:id:cross_hyou:20190926194544j:plain

県内総生産額は鳥取県が一番少ないです。11兆7960億35百万円です。

東京が一番多く102兆1295億28百万円です。

今回は以上です。

 

台風の発生数と上陸数のデータの分析4 - 平均気温、台風発生数、台風上陸数で並び替え、と思ったら元のデータがダメだった。

 

www.crosshyou.info

 の続きです。

今回は、前回作成したデータフレームの各変数を大きい順、小さい順に並び替えましょう。

 まずは、沖縄の平均気温の高い順に並び替えます。

 

Statistics: An Introduction Using R

Statistics: An Introduction Using R

 

 を参考にします。

f:id:cross_hyou:20190925192117j:plain

2019年の8月が一番高いですね! 次が1998年ですが。。。と思いましたがなんかおかしいですね。。。2019年のM7が4つもあります。。。

なんかもとのデータがダメですね。。

せっかく、いままで苦労したのに。。。

R言語で読み込んだときにちゃんと確認すべきでした。

方針を転換して、BornとLand、つまり台風の発生数と上陸数だけを分析の対象としましょう。

2つしか変数ないですから、データフレームよりもベクトルのほうが扱いやすいですね。hassnew, journewというベクトルがありますからこれで分析しましょう。

まずは、hassnew, journewをそれぞれカウントデータとして分析します。

上陸数を反応変数、発生数を説明変数にしましょう。

まずは、散布図を描きます。

f:id:cross_hyou:20190925193421j:plain

jitter関数をつかうことによってそれぞれのプロット点を微妙にずらしています。

両方とも0が圧倒的に多いです。発生数が多いほど上陸数は多いですね。

glm関数でfamily = poissonとして回帰分析します。

f:id:cross_hyou:20190925193742j:plain

hassnewの係数は0.4648でp値は2e-16よりも小さいです。hassnewが大きいほどjournewも大きいということですね。

回帰モデルの線を重ねてみます。

f:id:cross_hyou:20190925194155j:plain

f:id:cross_hyou:20190925194209j:plain

predict関数でモデルから推測される反応変数の値を取得することができます。

今回は以上です。

気温のデータがおかしかったのがショックです。分析する前にちゃんとチェックしないとダメですね。反省。

台風の発生数と上陸数のデータの分析3 - 気温のデータと台風の発生数、上陸数のデータを一つのデータフレームにまとめる。

 

www.crosshyou.info

 の続きです。前回はNAを0に変換しました。

今回は、台風の上陸回数、発生回数をベクトルにします。

 

まず、現状のデータフレームがどういうふうなのかhead関数で見てみます。

f:id:cross_hyou:20190921095145j:plain

HASSとJOURは1行目が1951年、2行目が1952年と横向きにデータが入っています。

kionは1行目が1951年1月、2行目が1951年2月と縦向きにデータが入っています。

kionのように縦向きにデータが入っているほうがいいので、HASSとJOURのデータを縦向きに変えましょう。

まず、HASSの中で必要なデータは2列目から13列目なので、その部分だけにします。

f:id:cross_hyou:20190921095737j:plain

このhassnewをas.matrix関数でマトリックスにします。

f:id:cross_hyou:20190921101606j:plain

t関数でhassnewの縦と横を転換します。

f:id:cross_hyou:20190921101803j:plain

これで、1列目が1951年、2列目が1952年、というようにデータがならびました。

as.vector関数でマトリックスがベクトルになります。

f:id:cross_hyou:20190921102116j:plain

できました。同じようにして、JOURもベクトルに変換します。

f:id:cross_hyou:20190921102518j:plain

これで気温のデータと発生数、上陸数のデータを合体させて一つのデータフレームにできます。kionの行数、hassnewとjournewのデータ数を確認しましょう。

f:id:cross_hyou:20190921103137j:plain

kionの行数、hassnewのデータ数、journewのデータ数、みんな816ですね。

それでは、data.frame関数でデータを一つにまとめましょう。

f:id:cross_hyou:20190921103717j:plain

ようやく、気温のデータと台風の発生数、上陸数のデータが一つのデータフレームにまとまりました。

今回は以上です。

 

台風の発生数と上陸数のデータの分析2 - R言語でNAを0に変換するには、is.na関数を使う。

 

www.crosshyou.info

 の続きです。前回はデータをR言語に読込んだところで終わりました。今回は、台風の発生数と上陸数のデータのNAを0に変換します。

R言語でNAを0に変換するには、is.na関数を使います。

vector <- c(1, 2, NA, 4, NA, 5)

というベクトル、vectorがあったら、

vector[is.na(vector)] <- 0

とすればいいはずです。

早速やってみます。

f:id:cross_hyou:20190919195800j:plain

うまくできました!

この作業をすべての列で実行すればいいですが、あと23回もありますから、関数を作って一括してNAを0に変換しましょう。

まず、ベクトルの中のNAを0に変換する関数を作ります。

f:id:cross_hyou:20190919200117j:plain

function関数で関数を作りました。それでは、hasseiのM2で確かめてみます。

f:id:cross_hyou:20190919200348j:plain

あれ?ダメですね。。どこが悪かったのかな?

f:id:cross_hyou:20190919200719j:plain

 

これでどうでしょうか?ベクトルを出すように1行追加しました。やってみます。

f:id:cross_hyou:20190919200918j:plain

やった!うまくいきました。最後にRのコンソールに結果を出力するようにしないといけないのですね。

こうして関数を作成したので、apply関数で全部の列に適用します。

f:id:cross_hyou:20190919201152j:plain

うまくNAを0に変換できました!zerotwo関数は全部の結果を出力してしまうので、すこし変更しましょう。

f:id:cross_hyou:20190919201423j:plain

最後の出力のときに、head関数ではじめの6つのデータしか出さないようにしました。

それでは、jourikuのデータフレームで実行しましょう。

f:id:cross_hyou:20190919201719j:plain

うまくできました。

こうして、NAを0に変換したので、summary関数で台風の発生数と上陸数のデータの平均値などを見てみましょう。

f:id:cross_hyou:20190919201912j:plain

あれれ!?  NAがあります。。。あ!zeroNAで処理した結果を改めてなにか別のデータフレームに保存しないといけなかったのですね。

もう一度やります。

f:id:cross_hyou:20190919202235j:plain

すばらしい! 台風の発生は年間では21から28回ですね。ひと月に6回というのが最大です。

上陸数も見てみます。

f:id:cross_hyou:20190919202530j:plain

とここまでやって、いま、間違いをおかしていることに気が付きました。zeroNA関数ははじめの6個のデータしか取得しない関数でした。

関数はzerotwo関数を使えばよかったですね。

f:id:cross_hyou:20190919203000j:plain

zerotwo関数で処理したデータフレーム、HASS(発生数)JOUR(上陸数)を作成しました。

発生数を見ると、年間で少ないときは14個、多いときは39個も台風が発生しています。

月間では10個が最多です。

上陸数も見てみましょう。

f:id:cross_hyou:20190919203207j:plain

上陸は年間で10回が最多です。月間では4回が最多ですね。

ということで今回は以上です。

 

台風の発生数と上陸数のデータの分析1 - 気象庁からデータを取得して、R言語に読込む

今回は、気象庁のホームページ

www.jma.go.jp


から、台風の発生数、上陸数のデータを取得して分析しようと思います。

f:id:cross_hyou:20190918193831j:plain

 

f:id:cross_hyou:20190918193843j:plain

 

f:id:cross_hyou:20190918193934j:plain

 

からデータをダウンロードしました。

台風の発生数、上陸数、那覇と東京の月別の平均気温を取得しました。

f:id:cross_hyou:20190918194026j:plain

気温はこんなデータファイルですね。

f:id:cross_hyou:20190918194118j:plain

上陸数はこんなファイルです。

R言語で処理しやすいようにファイルを加工してからread.csv関数で読込みます。

f:id:cross_hyou:20190918195200j:plain

こんな感じで台風発生回数のデータを読込み、hasseiという名前のデータフレーム保存しました。

f:id:cross_hyou:20190918195457j:plain

台風の上陸風はこんな感じです。

f:id:cross_hyou:20190918195714j:plain

気温データはこんな感じです。

まず、NAを0に変換しましょう。どうしたらいいのかな?

今回は以上です。次回にNAを0に変換しましょう。

 

 

読書記録 - 「ゲーム理論入門の入門」 鎌田雄一郎 著 岩波新書

 

ゲーム理論入門の入門 (岩波新書)

ゲーム理論入門の入門 (岩波新書)

 

 ゲーム理論入門の入門ということですが、難しかったです。

P37の図2-3のような図では、自分はよくわからなかったです。

P29のような利得表のほうがわかりやすいです。

最後のほうの完全ベイジアン均衡は理解できませんでした。

罰金ゲームの話が一番面白かったです。

 

日経の経済指標からデータを取得して分析4 - R言語でロジスティック回帰分析、そしてドル円だけが残った。

 

www.crosshyou.info

 の続きです。

今回は日経平均の騰落を回帰分析したいと思います。ロジスティック回帰分析です。

まずは、日経平均の騰落を表すベクトルを作成します。

f:id:cross_hyou:20190914154439j:plain

まず、vNikkの先頭にNAを追加したベクトル、xを作成しました。

これで、vNikk - x で前月比になります。

f:id:cross_hyou:20190914154856j:plain

先頭がNAなので削除します。

f:id:cross_hyou:20190914155016j:plain

できました。

説明変数のvCons, vIP, vYenも先頭を削除してデータの時点を合わせます。

f:id:cross_hyou:20190914155352j:plain

vToraをマイナスなら0, プラスな1に変換します。ifelse関数を使います

f:id:cross_hyou:20190914155937j:plain

平均値が0.58ですから、分析期間中の日経平均は、前月比プラスの月が58%あったということがわかります。

これで準備ができました。R言語でロジスティック回帰分析をするには、glm関数を使います。

 

Statistics: An Introduction Using R

Statistics: An Introduction Using R

 

 この本を参考にして分析します。

f:id:cross_hyou:20190914160426j:plain

vCons:vIP:vYenは削除してもよさそうですね。

f:id:cross_hyou:20190914160623j:plain

ロジスティック回帰分析のモデルをanova関数で比較するときは、test = "Chi"で比較するようです。P値が0.7766ですから、model2を採用します。

f:id:cross_hyou:20190914160832j:plain

vIP:vYenを削除してみます。

f:id:cross_hyou:20190914161018j:plain

p値が0.445と0.05より大きいですから、model3とmode2は有意な差は無いです。なので、model3を採用します。

f:id:cross_hyou:20190914161256j:plain

vCons:vYenの係数のp値は0.0937と0.05よりも大きいですから、これを削除してみます。

f:id:cross_hyou:20190914161513j:plain

p値が0.06238と0,05よりも大きいので、より単純なモデルのmodel4を採用します。

model4を調べます。

f:id:cross_hyou:20190914161712j:plain

vCons:vIPを削除してみましょう。

f:id:cross_hyou:20190914161900j:plain

p値が0.1792と0.05よりも大きいのでmodel5をさらに調べます。

f:id:cross_hyou:20190914162101j:plain

p値が一番大きいvConsを削除しましょう。

f:id:cross_hyou:20190914162243j:plain

p値が0.6404なので、model6を採用し、さらに調べます。

f:id:cross_hyou:20190914162415j:plain

I(vYen^2)を削除しましょう。

f:id:cross_hyou:20190914162607j:plain

p値が0.6854なのでmodel7を調べます。

f:id:cross_hyou:20190914162727j:plain

vIPの係数のp値が0.604と一番大きいので削除します。

f:id:cross_hyou:20190914162920j:plain

p値が0.6032なので、model8を採用して調査を続けます。

f:id:cross_hyou:20190914163057j:plain

I(vCons^2)を削除します。

f:id:cross_hyou:20190914163238j:plain

p値が0.2558なので、model9を採用します。model9を見てみましょう。

f:id:cross_hyou:20190914163401j:plain

I(vIP^2)を削除します。

f:id:cross_hyou:20190914163528j:plain

p値が0.2087と0.05より大きいですから、model10を採用します。model10を見てみましょう。

f:id:cross_hyou:20190914163703j:plain

あ!vYenの係数のp値が0.05よりも大きくなってしまいました!vYenも削除しましょう。

f:id:cross_hyou:20190914163905j:plain

p値が0.04581と0.05よりも小さいですから、vYenを削除するのはよくないということですね。

結局model10が最後のモデルです。

model10のvYenの係数の符号がマイナスですから、円安水準のほうが日経平均は上昇し易いということですね。

グラフで表現してみましょう。

f:id:cross_hyou:20190914164650j:plain

f:id:cross_hyou:20190914164705j:plain

青い線が回帰直線でですが、ダメですね。

株式相場の騰落を予測するのは難しいです。

今回は以上です。