www.crosshyou.info

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

都道府県別の生活保護被保護実世帯数データの分析1 - R言語でCSVファイルのデータを読み込む。大阪府が一番多い。福井県が一番少ない。

今回は、都道府県別の生活保護被保護実世帯数データを分析してみようと思います。

データは、政府統計の総合窓口e-Statから取得します。

f:id:cross_hyou:20200109191037p:plain

まず、47都道府県を選択し、

f:id:cross_hyou:20200109191110p:plain

総人口、可住地面積、県内総生産額と生活保護被保護実世帯数の4つのデータを選択します。

f:id:cross_hyou:20200109191203p:plain

CSVファイルはこんな感じです。9行目に変数名を追加しています。

これをR言語のread.csv関数で読込みます。

f:id:cross_hyou:20200109191836p:plain

str関数でデータフレームの構造を確認しました。うまく読込んでいます。

na.omit関数でNAのある行を削除しましょう。

f:id:cross_hyou:20200109192109p:plain

Prefのところの数を見ると10とありますから、10年間分のデータがあるようですね。

yearのファクタ水準を整理します。

f:id:cross_hyou:20200109192517p:plain

2006年度から2015年度までの10年間のデータがあるのですね。

このデータで私の興味のあることは、生活保護被保護実世帯数が人口、可住地面積、県内総生産額の3つの変数と関連しているかどうか、生活保護被保護実世帯数の増加(減少)がこれらの3つの変数と関連しているかどうか、の2点を確認したいと思います。

まずは、各変数の10年間の平均値を計算してみましょう。

f:id:cross_hyou:20200109193317p:plain

tapply関数で作成できます。生活保護被保護実世帯数をグラフにしてみましょう。

f:id:cross_hyou:20200109193553p:plain

f:id:cross_hyou:20200109193640p:plain

上方の外れ値がたくさんありますね。

小さい順に表示してみましょう。

f:id:cross_hyou:20200109193902p:plain

東京都よりも大阪府のほうが数が多いことがわかります。一番すくないのは福井県です。

今回は以上です。

 

日銀短観2019年12月調査のデータ分析7 - R言語で企業の規模と現状の景況感で先行きの景況感を分析。(ANCOVA)

 

www.crosshyou.info

の続きです。

今回は、企業の規模というカテゴリカル変数と現状の景況感という実数の変数を組み合わせて先行きの景況感を分析してみたいと思います。

ANCOVA(ANalysis of COVAriance)というものですね。

まずは、横軸に現状の景況感(Now)、縦軸に先行きの景況感(Next)の散布図を作成してみます。その際、企業の規模で色分けしてみます。

f:id:cross_hyou:20191228105454p:plain

f:id:cross_hyou:20200108195004p:plain

赤が大企業、緑が中堅企業、青が中小企業です。中小企業は大企業、中堅企業と比べるとNowもNextも水準が低いように見えます。

lm関数でANCOVAを実行します。

f:id:cross_hyou:20200108195620p:plain

NowとSizeの交互作用は有意ではないようです。交互作用が有意では無いということは、回帰線の傾きは企業規模で変らないということですね。削除してモデルを単純化します。

f:id:cross_hyou:20200108195925p:plain

model2のサマリを見てみます。

f:id:cross_hyou:20200108200232p:plain

中堅企業のEstimateのは-5.33355で、中小企業のEstimateは-4.54664です。差は0.8ぐらいしかありません。Standard Errorは1.5ぐらいありますから、中堅企業と中小企業とは有意な差はなさそうですね。新しいファクタを作って確認します。

f:id:cross_hyou:20200108200759p:plain

このnewSizeという新しいファクタでモデルを作成してみます。

f:id:cross_hyou:20200108201011p:plain

model2とmodel3では有意な違いは無いですね。

model3を確認します。

f:id:cross_hyou:20200108201142p:plain

p-valueが2.2e-16よりも小さいので有意なモデルです。先ほどのように、散布図を描いてみます。

f:id:cross_hyou:20200108201642p:plain

f:id:cross_hyou:20200108201625p:plain

この散布図にそれぞれの回帰直線を重ねてみます。

大企業の回帰直線は、切片が0.56625で傾きが0.81183です。

中堅中小の回帰直線は、傾きは同じ0.81183ですが、切片が0.56625 - 4.96150 = -4.39325になります。

f:id:cross_hyou:20200108202453p:plain

f:id:cross_hyou:20200108202544p:plain

こうなります。赤が大企業で、緑が中堅中小企業です。

最後にmodel3の残差プロットなどを見てみます。

f:id:cross_hyou:20200108203857p:plain

f:id:cross_hyou:20200108203908p:plain

今回は以上です。

 

読書記録 - 「オスマン帝国 繁栄と衰亡の600年史」小笠原弘幸 著 (中公新書)

 

 

オスマン帝国-繁栄と衰亡の600年史 (中公新書)

オスマン帝国-繁栄と衰亡の600年史 (中公新書)

  • 作者:小笠原 弘幸
  • 出版社/メーカー: 中央公論新社
  • 発売日: 2018/12/19
  • メディア: 新書
 

私の理解は間違っているかもしれないが、

オスマン帝国は600年もの間、ずっと男子の系統の一族が統治していた。その理由は兄弟を殺すことと、お嫁さんを奴隷から娶り、ハレムを作ることで親族内での後継者争いの種を事前に摘み取ることだった。

逆に言えばこうでもしないと600年の長きにわたり、ある一つの一族がずっと統治し続けることは不可能だと思った。

 

日銀短観2019年12月調査のデータ分析6 - 企業の規模(大企業/中堅企業/中小企業)でANOVA。企業規模によって景況感は違いがある。

 

www.crosshyou.info

 の続きです。

今回は、企業の規模、大企業、中堅企業、中小企業の3つの企業規模で短観の値に違いがあるのかどうか、ANOVA(ANalyis Of VAriance)をしてみます。

まずは、現状(Now)のそれぞれの平均値を見てみましょう。tapply関数とmean関数を使います。

f:id:cross_hyou:20191228094904p:plain

大企業の平均は9.92, 中堅企業の平均は9.61, 中小企業の平均は-0.14と中小企業は著しく低いですね。plot関数で箱ひげ図を描きます。

f:id:cross_hyou:20191228095158p:plain

f:id:cross_hyou:20191228095209p:plain

箱ひげ図を見ると、大企業も中堅企業も中小企業もけっこう上下の幅が大きいですね。平均値だけ見ると、絶対に中堅企業と中小企業では違いがあるな、と感じましたが、箱ひげ図を見た印象はそんなに違わないんじゃないか、と感じました。

実際はどうなんでしょうか?aov関数でANOVAができるそうです。やってみます。

f:id:cross_hyou:20191228095747p:plain

p-value が0.0373と0.05よりも小さいので、Sizeの違いはNowの違いに影響していますね。

R言語では、aov関数で簡単にANOVAができますが、これを手作業(といっても計算はRでやります)でやったらどうなるでしょうか?練習と思ってやってみます。

Michael J. CrawleyのStatistics An Introduction using Rを参考にしてやってみます。

 

Statistics: An Introduction Using R (English Edition)

Statistics: An Introduction Using R (English Edition)

  • 作者:Michael J. Crawley
  • 出版社/メーカー: Wiley
  • 発売日: 2014/09/23
  • メディア: Kindle版
 

 まずは、SSYというのを計算します。これは、個々の観測値から全体の平均値を引いて、それを2乗して、合計したものです。

f:id:cross_hyou:20191228100548p:plain

SSYは21200.81となりました。この値は、aovのsummary画面のどこにあるかというと、

f:id:cross_hyou:20191228100821p:plain

Sum Sqの1654と19547の合計値です。

次に、SSEを計算します。SSEはthe error of sum of squaresです。大企業、中堅企業、中小企業それぞれで、SSYのような計算をして、合計します。

f:id:cross_hyou:20191228101945p:plain

このSSEの値、19546.791はaovの結果の

f:id:cross_hyou:20191228102049p:plain

ここですね、Residualsの行のSum Sqの値が、SSEです。

そして、Sizeの行のSum Sqの値、1654がSSA(the treatment sum of suqares)です。

SSY = SSA + SSEという関係式が成り立ちます。

全体のバラツキ = ファクタ内のバラツキ + 残りのバラツキ

という関係かな。なので、SSA = SSY - SSEです。

f:id:cross_hyou:20191228102808p:plain

そして、このSSAの自由度(Df)は、2です。大企業、中堅企業、中小企業と3つのファクターがあるから1引いて2です。SSA(Sum Sq)の1654を自由度2で割った値が、Mean Sqの827.0になります。

f:id:cross_hyou:20191228103100p:plain

SSE(Residuals)の自由度はいくつでしょうか?全体の観測値の数は84です。SSYの自由度は83で、SSAの自由度は2ですから、83-2=81になります。

SSYの自由度 = SSAの自由度 + SSEの自由度

という関係ですね。

なので、SSE(Residuals)のMean Sqは、

f:id:cross_hyou:20191228103559p:plain

241.3183です。

aovの表でSSA(Size)のMean SqとSSE(Residuals)のMean Sqを確認しましょう。

f:id:cross_hyou:20191228103912p:plain

そして、SSA(Size)のMean SqがSSA(Size)のvarianceで、SSE(Residuals)のMean SqがSSE(Residuals)のvarianceですね。この二つのvarianceの比を計算します。

f:id:cross_hyou:20191228104258p:plain

この3.427がaovの結果画面のF valueの値です。

このF_ratioが統計的に有意かどうか、これをpf関数で計算します。

f:id:cross_hyou:20191228104547p:plain

はい、この0.037263がaovの結果の一番右の列、Pr(>F)の値と一致しましたね。

これが手計算でのANOVAです。結構面倒ですね。これをaov関数で一瞬でできてしまうのがRの凄いところですね。

NowはSizeによって違いがあることがわかりました。Next(先行き)はどうでしょうか?

f:id:cross_hyou:20191228105009p:plain

Next(先行き)も中小企業は低いですね。

箱ひげ図を描きます。

f:id:cross_hyou:20191228105221p:plain

f:id:cross_hyou:20191228105229p:plain

箱ひげ図はこんな感じでした。大企業のバラツキは小さいですね。

それでは、aov関数でANOVAを実行します。

f:id:cross_hyou:20191228105454p:plain

p値は0.00727です。Nowのときよりも小さなp値ですね。Next(先行き)でも企業規模によって景況感は違うことがわかりました。

今回は以上です。

 

日銀短観2019年12月調査のデータ分析5 - 現状と先行きの回帰分析。

 

www.crosshyou.info

 の続きです。

今回は、現状と先行きで回帰分析をしてみようと思います。

先行きを反応変数、現状を説明変数にしてみます。

まずは、plot関数で散布図を描いてみましょう。

f:id:cross_hyou:20191221101555p:plain

f:id:cross_hyou:20191221101605p:plain

右肩上がりの散布図です。現状(Now)と先行き(Next)には正の相関があることがわかります。cor.test関数で相関係数を見てみます。

f:id:cross_hyou:20191221101958p:plain

相関係数は0.9099091とかなり高い相関です。p-value < 2.2e-16です。95%信頼区間の相関係数は0.8640785から0.9407773となっています。現状(Now)と先行き(Next)には確固たる相関関係があります。

R言語で回帰分析をするには、lm関数を使います。

f:id:cross_hyou:20191221102538p:plain

lm関数で回帰分析モデルを作ってsummary関数でモデルを表示します。

p-value < 2.2e-16と0.05以下ですから、この回帰分析モデルは有意です。

切片とNowの係数のp値はともに0.05以下ですので、どちらも意味があります。

Adjusted R-squaredが0.8258です。これはこのモデルが先行き(Next)の値の82.58%を説明できている、ということです。

この回帰分析モデル式は、

先行き(Next) = -2.85264 + 0.83000 * 現状(Now)

という式です。現状(Now)の値に0.83を掛けて、2.85264を引いた値が先行き(Next)です。

ここからは、Rクックブック Paul Teetor著を参考にしていろいろやってみます。

 

Rクックブック

Rクックブック

  • 作者:Paul Teetor
  • 出版社/メーカー: オライリージャパン
  • 発売日: 2011/12/22
  • メディア: 大型本
 

 

plot関数でこのモデルの残差プロットを見てみます。

f:id:cross_hyou:20191221103343p:plain

f:id:cross_hyou:20191221103311p:plain

 

特に問題がありそうな残差プロットではないです。

confint関数で回帰式の係数の信頼区間がわかります。

f:id:cross_hyou:20191221104121p:plain

Nowから一番低いNextを導出する式は、

Next = -4.2677804 + 0.746879 * Now

です。

その反対に、Nowから一番高いNextを導出する式は、

Next = -1.437496 + 0.931127 * Now

です。この二つの式の直線とmodel1の直線を散布図に重ね合わせてみましょう。

f:id:cross_hyou:20191221105156p:plain

f:id:cross_hyou:20191221105401p:plain

confint関数はmatrixオブジェクトになっていますので、confint(model1)[ , 1]で一列目の値を指定できます。

influence.measures関数で回帰分析モデルに大きな影響を与えている観測値がわかります。

f:id:cross_hyou:20191221111110p:plain

一番右のアスタリスクが大きな影響を与えている、という印です。このデータの番号を一つ一つ確認すればいいのですが、もっと効率的にやります。

まず。このinfluenceというオブジェクトには、$is.infで影響力があればTRUE、なければFALSEのマトリックスがあります。

f:id:cross_hyou:20191221111436p:plain

このinfmtxというマトリックスを行で合計します。全部FALSEだったら0になってその行の観測値は影響が小さい、ということで一つでもTRUEがあれば1以上の値になってその観測値は影響が大きいとわかります。rowSums関数でマトリックスの行を合計できます。

f:id:cross_hyou:20191221111837p:plain

as.logical関数でTRUE/FALSEのベクトルに変換します。

f:id:cross_hyou:20191221112111p:plain

このidxでデータフレームからTRUEの行だけ表示すればいいですね。

f:id:cross_hyou:20191221112349p:plain

上の7つが回帰モデルに大きな影響を与えています。

散布図にしてみましょう。

f:id:cross_hyou:20191221112934p:plain

f:id:cross_hyou:20191221112948p:plain

今回は以上です。

 

日銀短観2019年12月調査のデータ分析4 - 現状の変化幅と先行きの変化幅に有意な違いはあるのか? 変化幅には有意な違いは無し。

 

www.crosshyou.info

 の続きです。

前回は現状(Now)と先行き(Next)の水準そのものに有意な違いがあるかどうかを調べました。その結果、現状と先行きには有意な違いがあり、先行きのほうが低いということがわかりました。

今回は、変化幅で同じように分析してみます。

まずは、summary関数で二つを比較してみます。

f:id:cross_hyou:20191219192751p:plain

NowChgの平均値は-2.988、NextChgの平均値は-3.893です。Nextのほうが悪いようですね。

hist関数でヒストグラムを見比べてみます。

f:id:cross_hyou:20191219193228p:plain

f:id:cross_hyou:20191219193239p:plain

NextChgのほうが尖っている分布ですね。

boxplot関数で箱ひげ図を見比べます。

f:id:cross_hyou:20191219193618p:plain

f:id:cross_hyou:20191219193648p:plain

前回はt.test関数で、paired = TRUEというのをつけてt検定を実行しましたが、今回は、NextChg - NowChgのベクトルを作って、これが0より小さいかどうかを検定しましょう。

まず、ベクトルを作成します。

f:id:cross_hyou:20191219194151p:plain

f:id:cross_hyou:20191219194133p:plain

なんとなくですが、0より小さい値のほうが多いように見えます。

それでは、t.test関数を実行します。

f:id:cross_hyou:20191219194410p:plain

p-value = 0.3782と0.05よりも大きいので、差が0より小さいとは言えないですね。結論はNowChgとNextChgの平均値には有意な違いは無い、ということです。

普通のt.test関数でpaired = TRUEでも確認してみます。

f:id:cross_hyou:20191219194659p:plain

p-value = 0.3782と同じ結果です。

Wilcoxon's Rank-Sum Testもやってみます。wilcox.test関数ですね。

f:id:cross_hyou:20191219195043p:plain

p-value = 0.2393と0.05よりも大きいです。NowChgとNextChgで違いがあるとは言えないですね。

今回は以上です。

日銀短観2019年12月調査のデータ分析3 - 現状と先行きに有意な違いはあるのか?先行きのほうが現状よりも有意に悪い。

 

www.crosshyou.info

 の続きです。

今回は現状と先行きに統計的に有意な違いがあるのかどうかを調べます。

まずは、summary関数でデータフレームを確認しましょう。

f:id:cross_hyou:20191218190456p:plain

今回もMichael J. CrawleyのStatistics An Introduction Using R Second Edition

 

Statistics: An Introduction Using R (English Edition)

Statistics: An Introduction Using R (English Edition)

  • 作者:Michael J. Crawley
  • 出版社/メーカー: Wiley
  • 発売日: 2014/09/23
  • メディア: Kindle版
 

 を参考にします。Chapter 6 Two Samplesを主に参考にします。

データフレームのNowが現状で、Nextが先行きです。

Nowの平均値は、6.119です。Nextの平均値は、2.226です。この違いは統計的に有意な違いなのでしょうか?

まず、NowとNextのvarianceが同じと言えるかどうかを見てみます。

f:id:cross_hyou:20191218191002p:plain

var関数でvarianceがわかります。

Nowのvarianceは、255.4314です。Nextのvarianceは212.5386です。

このvarianceの比、255.4314 / 212.5386をF_ratioという変数名で保存します。

f:id:cross_hyou:20191218191454p:plain

F_ratioは、1.201812です。var関数を使っても手計算でも同じ値です。

このF_ratioが大きい値だとNowとNextのvarianceが違うとわかります。

NowとNextのデータ個数は何個でしょうか?length関数で数えましょう。

f:id:cross_hyou:20191218191825p:plain

84個ですね。ということは、自由度は、84 - 1 = 83です。

F_ratioの分母も分子も自由度は83です。α = 0.05, 自由度83, 83のときのthe critical value of Fを求めます。qf関数です。

 

f:id:cross_hyou:20191218193508p:plain

the critical value of Fは1.542242です。F_ratioは1.201812とこれよりも小さいので、NowとNextのvarianceに違いがるとは言えないです。

var.test関数でこれが簡単にできます。

f:id:cross_hyou:20191218193759p:plain

p-value=0.4041と0.05よりも大きいです。つまり帰無仮説:二つのvarianceのratioは1である、を棄却できません。つまり、統計的に有意な違いは無いということです。

では、NowとNextの平均値に有意な違いがあるかどうか、t.test関数でt検定をします。

このケースでは、IndustryがNowとNextを両方回答していますから、対になったサンプルですね。なのでt.test関数にpaired = TRUEを加えます。

f:id:cross_hyou:20191218194347p:plain

p-value=6.704e-07と0.05よりも小さいので、NowとNextの平均値には有意な違いがある、といことです。NowよりもNextのほうが平均値は小さかったですね。(Nowの平均値は、6.119です。Nextの平均値は、2.226です。)、つまり日本企業全体として現状よりも先行きが厳しいとみているのでしょうね。

wilcox.test関数でWilcoxon Rank-Sum Testを実行します。

f:id:cross_hyou:20191218195349p:plain

p-value = 7.986e-07ですので、NowとNextでは有意に違うことがこちらの関数でも確認できました。

今回は以上です。