crosshyou

主にクロス表(分割表)分析をしようかなと思いはじめましたが、あまりクロス表の分析はできず。R言語の練習ブログになっています。

景気ウォッチャー調査地域別(現状)のデータ分析2 - R言語のtapply関数で年別、月別、地域別に集計する。

 

www.crosshyou.info

 の続きです。

今回は、年ごと、月ごと、地域ごとに集計をしようと思います。tapply関数を使えばいいです。

まずは年別の平均値です。

f:id:cross_hyou:20200118160229p:plain

barplot関数で棒グラフにしてみます。

f:id:cross_hyou:20200118160544p:plain

f:id:cross_hyou:20200118160600p:plain

abline関数で50の水準に赤い水平線を引きました。50以上だと景気が良くて、50未満だと景気が悪いことを意味します。2008年、リーマンショックの時はとても景気が悪かったですね。2018年、2019年と景気が悪くなってきています。

月別ではどうでしょうか?

f:id:cross_hyou:20200118160854p:plain

9月が一番大きい値のようです。それでも50以下ですが。。

barplot関数で棒グラフにしてみます。

f:id:cross_hyou:20200118161247p:plain

f:id:cross_hyou:20200118161302p:plain

こんどはabline関数で平均値の水準に水平線を引いてみました。3、4、5月が平均以下ですね。

地域別でも見てみましょう。

f:id:cross_hyou:20200118161525p:plain

sort関数を使って小さい順に表示しました。北関東が一番低く、沖縄が一番高いです。

barplot関数で棒グラフにします。

f:id:cross_hyou:20200118161911p:plain

f:id:cross_hyou:20200118161928p:plain

こんな感じです。東京都は平均以上なのに、北関東、関東、南関東は平均以下です。

tapply関数は2つのデータで集計することもできます。

f:id:cross_hyou:20200118163447p:plain

round関数を使って小数点以下2桁までにしました。

今回は以上です。

 

景気ウォッチャー調査地域別(現状)のデータ分析1 - データをR言語で読み取り、分析しやすいデータフレームにする。

今回は景気ウォッチャーのデータの分析をしてみようと思います。

内閣府のウェブからデータを取得しました。

f:id:cross_hyou:20200118132025p:plain

ここからエクセルをダウンロードして、今回は、地域別(現状)のデータを使います。

R言語に読込ませるように、少し加工しました。

f:id:cross_hyou:20200118132209p:plain

これをR言語で読込みます。

f:id:cross_hyou:20200118132843p:plain

こうなりました。

このような形式のデータはよくないですよね。本当は、

f:id:cross_hyou:20200118133146p:plain

こういう形式になってないといけないですよね。

どうやって変換するのかな。。。

str関数の結果を見ると、obsは216あります。216/12=18なので2002年から2019年までの18年間のデータですね。variablesが17で、そのうち二つはYearとMonthですから地域は15地域あります。ということは、Yearは一つの年で、12か月 x 15地域 = 180必要です。1年につき180で2002年から2019年まで18年ですから、データの長さは、180 x 18 = 3240です。rep関数をYEARという変数を作成します。

f:id:cross_hyou:20200118134804p:plain

これでいいですね。

つぎは、MONTHという変数をつくりましょう。1月が15地域、2月が15地域、、、11月が15地域、12月が15地域、これを18回繰り返します。

f:id:cross_hyou:20200118135341p:plain

これでいいですね。

REGIONは15の地域を3240 / 15 = 216回繰り返します。

f:id:cross_hyou:20200118135835p:plain

これでいいですね。

最後の景気ウォッチャーの値ですが、これは横方向にデータを並べればいいですね。

でも、Rは普通は縦方向にデータを読込みますので、縦横を変換する必要があります。

f:id:cross_hyou:20200118140309p:plain

いいですね。

これで、分析のためのデータフレームができます。

f:id:cross_hyou:20200118140638p:plain

できましたね。

この新しいデータフレーム、dfaが正しく作成されているか確認してみましょう。

もともとのdfで地域別の平均値を出してみます。apply関数とmean関数を使います。

f:id:cross_hyou:20200118140955p:plain

新しいデータフレーム、dfaから地域別の平均値を計算します。tapply関数とmean関数を使います。

f:id:cross_hyou:20200118141237p:plain

並びは違いますが、値は同じです。

これで、データ分析用のデータフレーム、dfaができました。

今回は以上です。

 

都道府県別の生活保護被保護実世帯数データの分析6 - R言語で世帯数の伸び率をダミー変数を加えて回帰分析

 

www.crosshyou.info

 の続きです。

今回は、世帯数の伸び率をR言語で回帰分析してみたいと思います。

まずは、どのようなデータか再確認します。summary関数を使います。

f:id:cross_hyou:20200118115559p:plain

最低で、1.143倍、最大で1.970倍、平均が1.459倍、中央値が1.517倍です。

hist関数でヒストグラムを描きます。

f:id:cross_hyou:20200118115915p:plain

f:id:cross_hyou:20200118115950p:plain

perGDP, perAreaとの散布図を見てみましょう。

f:id:cross_hyou:20200118120259p:plain

f:id:cross_hyou:20200118120309p:plain

左の散布図はperGDPがX軸ですが、右にポツンと点がありますね。これは東京の点ですね。東京ダミーという変数を用意してこれも回帰分析の変数に加えましょう。

右の散布図では右端の点は北海道ですから、北海道ダミーというのも加えましょう。

f:id:cross_hyou:20200118120851p:plain

これで回帰分析をします。lm関数を使います。

f:id:cross_hyou:20200118121111p:plain

p-valueが0.1314なのでモデル自体が有意ではないですね。生活保護の伸びは一人当りのGDPや一人当りの可住地面積とは関係なさそうですね。とりあえず、perGDP:perAreaを削除してもっと単純なモデルを作成します。update関数を使います。

f:id:cross_hyou:20200118121422p:plain

anova関数でmodel1とmodel2を比較しました。Pr(>F)が0.5175と0.05よりも大きいですから、model1とmodel2で有意な違いは無いです。model2を見てみます。

f:id:cross_hyou:20200118121625p:plain

p-valueは0.08639と0.05よりは大きいので有意なモデルではないです。でも、model1よりはp-valueは小さくなりました。perGDPを削除してみます。

f:id:cross_hyou:20200118121846p:plain

model2とmodel3では有意な違いは無いです。model3を見てみましょう。

f:id:cross_hyou:20200118122002p:plain

p-valueが0.04901と0.05以下になりました。有意なモデルになりました。TKがいらないようです。削除します。

f:id:cross_hyou:20200118122154p:plain

model3とmodel4で有意な違いはありません。model4を見てみます。

f:id:cross_hyou:20200118122325p:plain

p-valueは0.02167なので有意なモデルです。HKADはいらないようです。削除しましょう。結局、東京ダミーも北海道ダミーも必要なかったですね。

f:id:cross_hyou:20200118122519p:plain

model4とmodel5では有意な違いはありません。model5を見てみます。

f:id:cross_hyou:20200118122701p:plain

p-valueは0.00747ですから1%以下の水準で有意な統計モデルです。perAreaの係数がマイナスですので、perAreaが大きいほど、伸び率は小さいということです。

散布図と回帰直線を重ねてみます。

f:id:cross_hyou:20200118123028p:plain

f:id:cross_hyou:20200118123038p:plain

なんか、直線で回帰するより、曲線で回帰するほうがいいような気がします。

念のため確かめます。

f:id:cross_hyou:20200118123309p:plain

Pr(>F)が0.8127ですから、2乗項を追加しても有意な違いにはならないですね。

今回は以上です。

 

都道府県別の生活保護被保護実世帯数データの分析5 - R言語でGeneralized Additive Model

 

www.crosshyou.info

 の続きです。

今回は、gam関数でgeneralized additive modelというのをやってみたいと思います。

Michael J. Crawley の Statistics: An introduction using R

 

Statistics: An Introduction Using R

Statistics: An Introduction Using R

  • 作者:Michael J. Crawley
  • 出版社/メーカー: Wiley
  • 発売日: 2014/11/24
  • メディア: ペーパーバック
 

 を参考にしました。

f:id:cross_hyou:20200116190508p:plain

f:id:cross_hyou:20200116190526p:plain

summary関数でmodelを見てみます。

f:id:cross_hyou:20200116190802p:plain

こうなりましたが、何を見たらいいのか全然わかりません。

前回のlm関数で作成したmodel2と比較したいと思います。AICを比較します。

f:id:cross_hyou:20200116193430p:plain

AICはmodel2のほうが小さいですから、lm関数で分析した線形モデルのほうがいいです。

今回は以上です。

 

都道府県別の生活保護被保護実世帯数データの分析4 - R言語で人口当りのデータで回帰分析

 

www.crosshyou.info

 の続きです。

前回は生活保護被保護実世帯数そのものを反応変数にして回帰分析をしました。

今回は、人口のデータで割って、人口当りのデータに直して回帰分析してみましょう。

まずは、各変数を人口で割ります。

f:id:cross_hyou:20200115192037p:plain

人口は1000を掛けているので、1000人当りの世帯数です。

GDPは100万円単位なので、人口一人当りのGDP(百万円)になります。

面積は100を掛けているので、人口一人当りの面積(100ha)です。

それではsort関数で並び替えて表示してみます。

f:id:cross_hyou:20200115192414p:plain

大阪府が一番多いですね。富山県が一番少ないです。

 

人口一人当りのGDPはどうでしょうか?

f:id:cross_hyou:20200115192932p:plain
東京都が一番多く、奈良県が一番少ないです。

 

人口一人当りの面積(100ha)はどうでしょうか?

f:id:cross_hyou:20200115193335p:plain

北海道が一番広く、東京都が一番狭いです。

 

各変数間の散布図を描きます。plot関数ですね。

f:id:cross_hyou:20200115193649p:plain

f:id:cross_hyou:20200115193702p:plain

各変数間の相関はあまり無いようですね。cor関数でみてみます。

f:id:cross_hyou:20200115193849p:plain

GDPとAreaが負の相関ですね。人口密度が高いほうが、一人当りのGDPは大きいですね。

それでは、lm関数で回帰分析をしてみます。

f:id:cross_hyou:20200115194303p:plain

p-valueは6.851と0.05よりも低いので有意なモデルです。I(perGDP^2)はいらないようです。削除してモデルを単純化します。

f:id:cross_hyou:20200115194512p:plain

Pr(>F)の値が0.35と0.05よりも大きいので、modelとmodel2で有意な違いはありません。model2を見てみましょう。

f:id:cross_hyou:20200115194704p:plain

model2の式は、

保護世帯数 = 12.24 + 0.78 * perGDP + 0.83 * perArea + 0.019 * perArea^2 - 0.45 * perGDP * perArea

という式になります。

どの都道府県がモデルから推測される世帯数と乖離しているか見てみましょう。

f:id:cross_hyou:20200115195415p:plain

埼玉県が一番実際の世帯数が少なく、大阪府が一番多いです。

今回は以上です。

都道府県別の生活保護被保護実世帯数データの分析3- R言語で回帰分析

 

www.crosshyou.info

 今回はR言語のlm関数で回帰分析をしてみます。

反応変数は、生活保護被保護実世帯数(avgHogo)で、説明変数は人口(avgPop), 可住地面積(acgArea), 県内総生産額(avgGDP)です。

まずは、相関マトリックスを見てみます。

f:id:cross_hyou:20200111153212p:plain

人口(avgPop)と県内総生産(avgGDP)の相関係数は0.9177とかなり相関が強いので、まずはavgPopとavgAreaを説明変数にしてみます。

f:id:cross_hyou:20200111153552p:plain

p-valueは8.34e-16と0.05よりも小さいので統計的に有意なモデルです。

avgPop:avgAreaは有意ではないようなので、削除します。

f:id:cross_hyou:20200111153826p:plain

avgPop:avgAreaを削除した、model2で問題ないです。確認します。

f:id:cross_hyou:20200111154018p:plain

avgAreaを削除してもよさそうですね。やってみます。

f:id:cross_hyou:20200111154155p:plain

model3でいいようです。確認します。

f:id:cross_hyou:20200111154302p:plain

model3のp-valueは2.2e-16よりも小さいので有意ですね。

散布図と回帰直線を描いてみます。

f:id:cross_hyou:20200111154556p:plain

f:id:cross_hyou:20200111154607p:plain

こんどは、説明変数をavgPopからavgGDPに変えてやってみます。

f:id:cross_hyou:20200111154808p:plain

avgGDP:avgAreaは必要ないようですので削除しましょう。

f:id:cross_hyou:20200111155019p:plain

nmodel2を見てみましょう

f:id:cross_hyou:20200111155211p:plain

avgAreaの係数のp値が0.0166なので有意ですね。nmodel2自体のp-valueは3.87e-14と0.05以下で有意です。

説明変数がavgGDP, avgAreaと二つあるので、coplot関数で散布図を描いてみましょう。

f:id:cross_hyou:20200111155918p:plain

f:id:cross_hyou:20200111155930p:plain

最後は、説明変数にavgPop, avgArea, avgGDPの三つを入れたモデルを調べてみましょう。

f:id:cross_hyou:20200111160231p:plain

avgPop:avgArea:avgGDPを削除してみましょう。update関数を使います。

f:id:cross_hyou:20200111160505p:plain

fmodel2を見てみます。

f:id:cross_hyou:20200111160658p:plain

avgPop:avgGDPは削除してかまわないようです。

f:id:cross_hyou:20200111160852p:plain

fmdel3を見てみます。

f:id:cross_hyou:20200111161007p:plain

avgPop:avgArea, avgArea:avgGDPの項とavgGDPの項が有意ですね。

model3, nmodel2, fmodel3の中でどれが一番いいのかな?残差プロットを見比べてみましょう。

f:id:cross_hyou:20200111161615p:plain

f:id:cross_hyou:20200111161626p:plain

どうなんでしょうね。。正直よくわかりません。

AIC関数でAIC(Akaike's Information Criterion)を算出してみました。

f:id:cross_hyou:20200111161819p:plain

AICが一番小さいのはfmodel3ですね。

もう一度fmodel3を見てみます。

f:id:cross_hyou:20200111162427p:plain

生活保護被保護実世帯数は、人口、可住地面積、県内総生産の3つの変数と関連があることがわかりました。

今回は以上です。

 

都道府県別の生活保護被保護実世帯数データの分析2 - R言語で世帯数の伸び率を計算。全体で10年間で1.5倍以上の伸び。埼玉県が一番の伸び率。

 

www.crosshyou.info

 の続きです。

今回は、生活保護被保護世帯数の伸び率を計算してみようと思います。

まずは、全国合計の値を計算しましょう。

2006年度の世帯数の合計を計算します。

f:id:cross_hyou:20200111103551p:plain

107万8368世帯です。

2015年度も同じようにします。

f:id:cross_hyou:20200111103848p:plain

162万9754世帯です。1.5倍以上ですよ。。。私の想像以上に増えていました。

ちゃんと計算します。

f:id:cross_hyou:20200111104202p:plain

1.511315倍です。

それでは、都道府県別の伸び率を計算します。

まず、tapply関数を利用して、2006年度の都道府県別の世帯数のデータを作成します。

f:id:cross_hyou:20200111105415p:plain
idx <- df$year == "2006年度"

として2006年度ならTrue, そうでないならFalseの論理ベクトルを作っています。Hogo2006という名前のデータを作成しました。

同じように、2015年度のデータも作成します。

f:id:cross_hyou:20200111105541p:plain

Hogo2015というデータを作成しました。

それでは、伸び率を計算します。

f:id:cross_hyou:20200111105813p:plain

一番伸び率の小さいのは山口県で1.143倍です。すべての都道府県で生活保護を受けている世帯の数は増えているのですね。一番伸び率の大きいのは埼玉県で1.969倍です。約2倍も増えています。10年で2倍も増えるって凄いですね。

グラフにして伸び率の分布を見てみます。

f:id:cross_hyou:20200111110221p:plain

f:id:cross_hyou:20200111110237p:plain

ヒストグラムの形状を見ると、山のピークが二つあることがわかります。

最後に次回以降の分析のために、前回計算した都道県別の平均値と今回計算した伸び率のデータを一つのデータフレームにまとめましょう。

f:id:cross_hyou:20200111111040p:plain

今回は以上です。