crosshyou

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

都道府県別の老人福祉費と児童福祉費の分析9 - R言語で、老人福祉費 / 児童福祉費 という比率を計算。沖縄県以外はどこも老人福祉費のほうが多い。

 

www.crosshyou.info

 の続きです。

今回は、老人福祉費 / 児童福祉費 という比率を計算してみます。

まず、老人福祉費と児童福祉が同じ都道府県の順番で並んでいるか確認します。

f:id:cross_hyou:20191019104251j:plain

== で同じかどうかをテストしました。すべてTRUEなので同じですね。sumでTRUEの数を数えると、47ですから47都道府県すべて同じ順番です。

それでは比率を計算します。

f:id:cross_hyou:20191019104445j:plain

沖縄県は比率が1.0以下ですから老人福祉費よりも児童福祉費のほうが多いのですね。その他の都道府県はすべて比率は1.0以上です。新潟県が3.27で一番高い比率です。

summary関数で平均値や中央値を調べます。

f:id:cross_hyou:20191019104826j:plain

中央値は2.1780、平均値は2.1671です。老人福祉費が児童福祉費の2倍強というのが平均ですね。

3つのグラフを描いてみます。

f:id:cross_hyou:20191019105059j:plain

f:id:cross_hyou:20191019105109j:plain

山型の分布です。

それでは、この老人福祉費 / 児童福祉費, OldChildをavgPop, avgArea, avgGDPの3つの変数で重回帰分析をしてみましょう。

まずは、一番複雑なモデルからスタートします。

f:id:cross_hyou:20191019105553j:plain

p-valueが0.1715と0.05よりも大きいですから、このモデルは有意ではないということです。老人福祉費 / 児童福祉費 は人口、面積、県内総生産では説明がつかないということですね。一応、step関数でモデルを単純化してみます。

f:id:cross_hyou:20191019110138j:plain

以下省略して、summaryを表示します。

f:id:cross_hyou:20191019110209j:plain

avgPop:avgGDPのp値は0.0866ですから不要かもですね。確認します。

f:id:cross_hyou:20191019110458j:plain

Pr(>F)の値が0.08664と0.05よりも大きいですから、model3を採用します。

f:id:cross_hyou:20191019110639j:plain

p-valueは0.04536と0.05より小さいので有意なモデルです。avgPopはいらなさそうです。

f:id:cross_hyou:20191019110827j:plain

model3とmodel4に有意な差は無いです。model4を採用します。

f:id:cross_hyou:20191019110957j:plain

p-valueは0.03841と0.05より小さいので有意なモデルです。avgGDPもいらなさそうですね。

f:id:cross_hyou:20191019111149j:plain

p値は0.4378と0.05よりも大きいので、model5を採用します。

f:id:cross_hyou:20191019111307j:plain

p-valueは0.01448と0.05よりも小さいので有意なモデルです。このmodel5はavgAreaが変数の単回帰モデルですね。I(avgArea^2)を追加してみます。

f:id:cross_hyou:20191019111544j:plain

p値が0.1762と0.05よりも大きいので、I(avgArea^2)を追加しても意味ないということですね。model5が最終的なモデルです。

散布図と回帰直線を描きます。

f:id:cross_hyou:20191019111913j:plain

f:id:cross_hyou:20191019111923j:plain

散布図の一番右のプロット、北海道が大きく影響しているようですね。

北海道を除外したら結果はまた変わるかもしれません。もしくは対数値の人口、面積、県内総生産で回帰分析したら結果はかわるかもしれないですね。

今回は以上です。

 

都道府県別の老人福祉費と児童福祉費の分析8 - R言語で一人当り児童福祉費を重回帰分析。東京都を含めるか除外するかでモデルが違ってくる。

 

www.crosshyou.info

 の続きです。

今回は一人当りの児童福祉費(ChildpM)を重回帰分析します。

説明変数は対数をとった面積(logArea)と対数をとった県内総生産(logGDP)です。

まずは、一番複雑なモデルから。

f:id:cross_hyou:20191017191555j:plain

logArea:logGDPの交差項は不要のようです。削除します。

f:id:cross_hyou:20191017191759j:plain

p値は0.9481です。model1もmodel2も有意な違いは無いです。より単純なmodel2を調べます。

f:id:cross_hyou:20191017192007j:plain

I(logArea^2)は不要ですね。削除します。

f:id:cross_hyou:20191017192206j:plain

model2とmodel3では有意な違いはありません。model3をさらに調べます。

f:id:cross_hyou:20191017192335j:plain

logAreaは不要ですね。削除します。前回の一人当り老人福祉費ではlogAreaは有意な変数でしたが、一人当りの児童福祉費では都道府県の面積は関係無いことがわかりました。

f:id:cross_hyou:20191017192721j:plain

p値が0.9822です、model3とmodel4で有意な差はありません。

model4をみてみましょう。

f:id:cross_hyou:20191017192941j:plain

これで、切片、logGDP、I(logGDP^2)の3つとも有意なモデルです。

logGDPをいろいろな値にして予測値を算出してみます。

f:id:cross_hyou:20191017193718j:plain

f:id:cross_hyou:20191017193732j:plain

こうなりました。散布図の一番右端の点、東京都ですがこれが影響して回帰曲線が放物線状になっていますね。東京都を除外すれば、単純に直線になるのではないでしょうか?やってみましょう。

f:id:cross_hyou:20191017194306j:plain

まずは、上のようにして東京都を除外したベクトルを作成しました。

後は同じ手順です。まずは一番複雑なモデルからやりましょう。

f:id:cross_hyou:20191017194545j:plain

exTArea:exTGDPはいらないです。削除します。

f:id:cross_hyou:20191017194743j:plain

exTm2を採用します。

f:id:cross_hyou:20191017194859j:plain

I(exTArea^2)は不要ですね。

f:id:cross_hyou:20191017195050j:plain

exTm3を採用します。

f:id:cross_hyou:20191017195215j:plain

exTAreaは必要ないようです。

f:id:cross_hyou:20191017195355j:plain

exTm4を採用します。

f:id:cross_hyou:20191017195522j:plain

これで、切片、exTGDP、I(exTGDP^2)の3つが有意になりました。2乗項も残りましたね。散布図と回帰曲線を描きましょう。

f:id:cross_hyou:20191017200011j:plain

f:id:cross_hyou:20191017200022j:plain

最後に東京都有りの散布図と回帰曲線、東京都無しの散布図と回帰曲線を並べて描いてみます。

f:id:cross_hyou:20191017200501j:plain

f:id:cross_hyou:20191017200512j:plain

東京都を入れないほうが回帰曲線のあてはまりがいいようです。

実際に東京都有りのmodel4のAdjusted R-squaredは0.5681

東京都無しのexTm4のAdjusted R-squaredは0.6873と東京都無しのほうが高いです。

今回は以上です。

 

都道府県別の老人福祉費と児童福祉費の分析7 - R言語で重回帰分析。一人当りの老人福祉費は面積が大きい県ほど多い。GDPの大きい県ほど少ない。

 

www.crosshyou.info

 の続きです。

今回は一人当りの老人福祉費を面積(対数をとったもの)と県内総生産(対数をとったもの)の二つの変数で回帰分析してみたいと思います。

まずは、それぞれの変数との散布図を描いてみます。

f:id:cross_hyou:20191016190159j:plain

f:id:cross_hyou:20191016190211j:plain

面積は関係なさそうですが、GDPは関係ありそうですね。

lm関数で重回帰分析のモデルを作り、調べていきます。

まずは、交差項も2乗項も入ったmaximumモデルからスタートします。

f:id:cross_hyou:20191016190613j:plain

logArea:logGDPの交差項のp値は0.2071ですので削除しましょう。

f:id:cross_hyou:20191016190818j:plain

anova関数でmodel1とmodel2を比較しました。p値が0.2071ですのでmodel2でいいですね。中身をみてみましょう。

f:id:cross_hyou:20191016191037j:plain


I(logArea^2)は必要なさそうですね。削除してみましょう。

f:id:cross_hyou:20191016191223j:plain

p値が0.612ですから、model2とmodel3は有意な差は無いです。従いまして、より単純なmodel3を採用します。summary関数で見てみます。

f:id:cross_hyou:20191016191442j:plain

あ、これですべての変数が有意になりました。モデル全体のp-valueは2.482e-10ですから有意なモデルです。Adjusted R-squaredは0.643ですから、都道府県別の一人当り老人福祉費は、面積と県内総生産、県内総生産の2乗で64.3%が説明できると言えます。

logAreaの係数は5.790とプラスですから、面積が大きい県ほど、一人当り老人福祉費は多いということですね。logGDPに関しては、logGDPの係数は-147, 2乗の項の係数は9.866ですから、

9.866logGDP^2 - 147.078logGDP = 9.866logGDP(logGDP - 14.907569)

と変換できますので、14.907569の半分の7.5ぐらいまでは、logGDPが0kら7.5ぐらいまではだんだんと減少して、それからはだんだんと増加していく、ということですね。

Y =  9.866 * X^2 - 147.078のグラフを描いてみましょう。

f:id:cross_hyou:20191016193304j:plain

f:id:cross_hyou:20191016193342j:plain

logGDPは6.25から8.00の範囲ですからこの範囲で考えると、6.25から7.5ぐらいまではlogGDPが大きいほど、一人当り老人福祉費は少なく、7.5から8.0の間ではlogGDPが大きいほど一人当り老人福祉費も大きいということですね。

model3の残差プロットを描きましょう。

f:id:cross_hyou:20191016193821j:plain

f:id:cross_hyou:20191016193836j:plain

今回は以上です。

 

日銀の短観データの分析6 - R言語で棒グラフと1標準誤差、信頼区間を表示する。

 

www.crosshyou.info

 の続きです。

今回は棒グラフと信頼区間を表示してみたいと思います。

 

Statistics: An Introduction Using R

Statistics: An Introduction Using R

 

 を参考図書にしてやってみます。

まず、Sector別の短観の平均値をtapply関数とmean関数で計算します。

f:id:cross_hyou:20191014142151j:plain

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

f:id:cross_hyou:20191014142749j:plain

f:id:cross_hyou:20191014142803j:plain

次に信頼区間を表示する関数を作成します。これがMichael J. Crawleyの本からの関数です。

f:id:cross_hyou:20191014143403j:plain

上の関数でyが棒グラフの部分です。これは今回はsectorbetsuです。zが信頼区間の大きさです。今回は、1standard errorを使います。

aov関数とsummary関数で1standard errorを確認します。

f:id:cross_hyou:20191014143953j:plain

このANOVA表から127がpooled error varianceだとわかります。

こんどは、Sectorがそれぞれ何個あるかをtable関数で確認します。

f:id:cross_hyou:20191014144427j:plain

製造業は228、全産業は12、非製造業は156です。なので、それぞれの1standard errorは、

f:id:cross_hyou:20191014144806j:plain

この3つの値が1standard errorになります。

これで、1標準誤差の表示ができます。

f:id:cross_hyou:20191014145244j:plain

f:id:cross_hyou:20191014145255j:plain

あれれ?三つの標準誤差の線がそれぞれの棒に表示されてしまいまたね。。

あの関数はすべての棒グラフで同じデータ数でないとだめなんですね。

あの関数を改良します。

f:id:cross_hyou:20191014150720j:plain

zもz[i]としないといけなかったですね。

これやってみましょう。

f:id:cross_hyou:20191014151055j:plain

f:id:cross_hyou:20191014151029j:plain

できました!

ちなみに関数をつかわないでarrows関数でひとつひとつ追加するのはこちらです。

f:id:cross_hyou:20191014151413j:plain

f:id:cross_hyou:20191014151427j:plain

1標準誤差をグラフに表示できましたので、今度は95%の信頼区間を表示してみましょう。

まずは、製造業、全産業、非製造業のそれぞれの信頼区間を計算します。

f:id:cross_hyou:20191014152027j:plain

このciをerror_bars2関数に使います。

f:id:cross_hyou:20191014152729j:plain

f:id:cross_hyou:20191014152748j:plain

こうしてグラフで95%信頼区間を表示すると、製造業と全産業、全産業と非製造業には違いは無いかもしれませんが、製造業と非製造業には明らかな違いがあることがわかります。
今回は以上です。

 

日銀の短観データの分析5 - R言語でANOVA。繊維や紙・パルプは景気が悪く、対事業所サービスや通信は景気がいい。

 

www.crosshyou.info

 の続きです。

今回の説明変数はIndus, 業種です。summary関数でどういう業種があるか見てみます。

f:id:cross_hyou:20191010154900j:plain

33種類の業種があります。今回は前回までとは違ったアプローチでANOVAをやってみたいと思います。

いつものように参考図書は、

 

Statistics: An Introduction Using R

Statistics: An Introduction Using R

 

 です。

aov関数でANOVAのモデルを作成し、summary.lm関数で表示します。

f:id:cross_hyou:20191010160106j:plain

f:id:cross_hyou:20191010160119j:plain

各業種のp値を見ると、電気・ガスが0.875278と一番大きいです。電気・ガスは、はん用機械と一緒にしても大丈夫そうです。

f:id:cross_hyou:20191010161152j:plain

p値が0.8753なのでmodel2でもいいですね。summary.lm関数で確認します。

f:id:cross_hyou:20191010161452j:plain

f:id:cross_hyou:20191010161508j:plain

非製造業が0.6878のp値で一番高いので、これを、はん用機械と統合してmodel3を作りましょう。

f:id:cross_hyou:20191010162138j:plain

p値が0.6878と0.05以上ですから、model2とmodel3に有意な違いはありません。なのでシンプルなほう、model3を採用します。model3を見てみましょう。

f:id:cross_hyou:20191010162644j:plain

f:id:cross_hyou:20191010162657j:plain

運輸・郵便のp値が一番大きいようです。これをはん用機械に統合してmodel4を作成します。

f:id:cross_hyou:20191010163219j:plain

p値が0.4799と0.05よりも大きいので、model3とmodel4に有意な違いはないです。なので、model4を採用してさらに統合できる業種は統合します。

f:id:cross_hyou:20191010163710j:plain

f:id:cross_hyou:20191010163734j:plain

化学のp値が一番大きいので化学をはん用機械に統合したmodel5を作成します。

f:id:cross_hyou:20191010164151j:plain

p値は0.3896と0.05より大きいので、model4とmodel5で有意な違いはあありません。従いまして、よりシンプルなmodel5を採用します。

f:id:cross_hyou:20191010164713j:plain

f:id:cross_hyou:20191010164727j:plain

業務用機械のp値が一番大きいので、これをはん用機械に統合します。

f:id:cross_hyou:20191010165147j:plain

p値は0.3351なのでmodel5とmodel6のふたつに有意な違いはありません。model6を調べます。

f:id:cross_hyou:20191010165508j:plain

f:id:cross_hyou:20191010165522j:plain

窯業・土石製品をはん用機械に統合します。

f:id:cross_hyou:20191010165914j:plain

p値が0.3665なので、model6とmodel7では有意な違いはありません。model7を採用します。

f:id:cross_hyou:20191010170217j:plain

f:id:cross_hyou:20191010170231j:plain

全産業をはん用機械に統合します。

f:id:cross_hyou:20191010170852j:plain

model8を採用します。

f:id:cross_hyou:20191010171107j:plain

f:id:cross_hyou:20191010171120j:plain

生産用機械をはん用機械に統合します。

f:id:cross_hyou:20191010171503j:plain

model9を採用します。

f:id:cross_hyou:20191010171852j:plain

f:id:cross_hyou:20191010171909j:plain

木材・木製品をはん用機械に統合します。

f:id:cross_hyou:20191010172215j:plain

model10を採用します。

f:id:cross_hyou:20191010172444j:plain

f:id:cross_hyou:20191010172501j:plain

卸売をはん用機械に統合します。

f:id:cross_hyou:20191010173041j:plain

model11を採用します。

f:id:cross_hyou:20191010173401j:plain

f:id:cross_hyou:20191010173415j:plain

とうとう、各業種が全部有意になりました。これでモデルの単純化は終わりです。

切片の値は8.6742です。これははん用機械と統合された業種の平均値が8.6742ということですね。物品賃貸は14.7424となっています。これは、8.6742 + 14.7424 = 23.416が物品賃貸の平均値ということですね。tapply関数でたしかめます。

f:id:cross_hyou:20191010173920j:plain

たしかにそうなっていますね。

最後にmodel11のグラフを描いておわります。

f:id:cross_hyou:20191010174232j:plain

f:id:cross_hyou:20191010174244j:plain

 

日銀の短観データの分析4- R言語でANOVA。中小企業は景況感は悪い。

 

www.crosshyou.info

 の続きです。

今回は、大企業、中堅企業、中小企業という企業規模の違い、Scaleをexplanatory variableにしてANOVAをしてみます。企業規模の違いで短観の数値に違いはあるでしょうか?

まずは、グラフで様子を確認しましょう。

f:id:cross_hyou:20191009231137j:plain

黒い点が大企業で黒い水平線が大企業の平均値。

赤い点が中堅企業で赤い水平線が中堅企業の平均値。

緑の点が中小企業で緑の水平線が中小企業の平均値です。

企業規模が小さいほど短観の値は悪くなっています。

tapply関数でそれぞれの平均値を計算します。

f:id:cross_hyou:20191009231729j:plain

大企業が11.07, 中堅企業は8.82, 中小企業は1.42と中小企業は極端に低いですね。

今回はtapplyの結果をkiboというベクトルに保存しました。なぜかというと、あとて、kibo[1]で大企業の平均値、kibo[2]で中堅企業の平均値、kibo[3]で中小企業の平均値として使うためです。

f:id:cross_hyou:20191009232047j:plain

ついでに、df1$Scale == "大企業", df1$Scale == "中堅企業”, df1$Scale == "中小企業"もベクトルにしておいてあとで、SSEの計算で使います。

f:id:cross_hyou:20191009232425j:plain

SSY, the sum of the squares of the differences between the y values and the allovermean を計算します。

f:id:cross_hyou:20191009232718j:plain

SSY, 毎回計算していますが、本当は1回計算すればいいです。でもこれは練習ブログなので毎回計算して、毎回同じ68298.76を得ています。

SSE, the error sum of squaresを計算します。ここでさきほど作った、kibo, idxL, idxM, idxSが活躍します。

f:id:cross_hyou:20191009233322j:plain

どうでしょう。前回はそれぞれのファクタのSSEを計算して最後に合計しましたが、今回はいちどにSSEを計算しました。わかりやすい式だと思います。SSEは61578.27です。
SSAはSSYからSSEを引いた値です。

f:id:cross_hyou:20191009233613j:plain

SSAは6720.49です。SSAのDegrees of freedom, 自由度は2です。なのでSSAのMean squareは、

f:id:cross_hyou:20191009233802j:plain

3360.245です。

SSEの自由度は、SSYの自由度が395で、SSAの自由度が2ですから、395 - 2 = 393です。SSEのMean squareは、

f:id:cross_hyou:20191009234039j:plain

156.6877です。なので、F ratioは3360.245 / 156.6877です。

f:id:cross_hyou:20191009234210j:plain

21.44549です。

これでANOVA表の構成要素は全部計算しました。

まとめましょう。

f:id:cross_hyou:20191009234846j:plain

F ratioの21.4という値が偶然に起こりえる値なのか、それとも企業規模に違いが無いと怒らない値なのか、qf関数でチェックします。

f:id:cross_hyou:20191009235128j:plain

3.02よりも21.4は大きいですから企業規模の違いは短観の値に影響があるといえます。

pf関数でp値も求めます。

f:id:cross_hyou:20191009235348j:plain

1.446713e-09です。0.05よりも小さい値です。

aov関数とsummary関数でいままでの答え合わせをします。

f:id:cross_hyou:20191009235610j:plain

p値もpf関数で計算した値と同じだし、ANOVA表ともSum of squaresやMean square, F ratioが一緒です。

今回は以上です。

企業規模は短観の数値と関係し、規模が小さくなるほど短観の値は低くなっています。

 

日銀の短観データの分析3 - R言語でANOVA。製造業よりも非製造業のほうが景気はいいようだ。

 

www.crosshyou.info

 の続きです。今回もANOVAをします。今度は、Sector(製造業、非製造業、全産業)を説明変数にして、Value、短観の値を反応変数にします。

またグラフでSectorによって短観の値に違いがあるか見てみましょう。

f:id:cross_hyou:20191009160101j:plain

f:id:cross_hyou:20191009160113j:plain

黒い点が製造業で黒い水平線が製造業の短観の平均値です。

赤い点が全産業で赤い水平線が全産業の短観の平均値です。

緑の点が非製造業で緑の水平線が非製造業の平均値です。業種による違いは大きそうですね。

tapply関数でそれぞれの業種の平均値を調べましょう。

f:id:cross_hyou:20191009160731j:plain

製造業の平均が1.39, 全産業の平均が7.57, 非製造業の平均が15.4です。業種による違いはありそうですね。

SSY, sum of the squares of the differences between y values and the overall meanを計算します。

f:id:cross_hyou:20191009161235j:plain

SSYは68298です。

SSE, the error sum of squaresを計算します。

f:id:cross_hyou:20191009162205j:plain

SSEは50073です。

SSAはSSY - SSEです。

f:id:cross_hyou:20191009162440j:plain

ここまでの結果をもとに、ANOVA表を作成します。

f:id:cross_hyou:20191009164718j:plain

D3のセル、ErrorのMean squareを計算します。

f:id:cross_hyou:20191009164932j:plain

f:id:cross_hyou:20191009165059j:plain

E2のセル、F ratioを計算します。

f:id:cross_hyou:20191009165201j:plain

f:id:cross_hyou:20191009165901j:plain

このF ratioの143が意味のある大きさかどうかが問題です。qf関数で調べます。

f:id:cross_hyou:20191009170033j:plain

143は3.865よりも大きいので、Sectorによって短観の結果に違いがある、ということですね。

aov関数とsummary関数で確認します。

f:id:cross_hyou:20191009170339j:plain

あ!Sectorは製造業、全産業、非製造業でしたから、Degrees of freedomは1ではなくて2ですね。さっきのANOVA表を訂正します。

f:id:cross_hyou:20191009170542j:plain

Mean square, F ratioを計算しなおします。

f:id:cross_hyou:20191009170808j:plain

f:id:cross_hyou:20191009170902j:plain

あらためて、aov関数とsummary関数の結果を提示します。

f:id:cross_hyou:20191009171006j:plain

p値が2e-16よりも小さいので、Sectorによって短観の値に違いはあります。

製造業よりも非製造業のほうが景況感はいいようですね。

今回は以上です。