crosshyou

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

都道府県別の通院者率のデータの分析5 - ダミー変数(year)の交差項を入れて回帰分析をする。

f:id:cross_hyou:20210923082217j:plain

Photo by Aaron Burden on Unsplash 

www.crosshyou.info

の続きです。

前回はhosp: 人口1000人当たりの通院者率をoldr: 65歳以上人口割合(%)で回帰分析しました。oldrが1ポイント高くなると、約8人通院者率が増えることがわかりました。

今回は回帰式にyear: 2007年と2019年のファクターを入れてみて、2007年と2019年では違いがあるかどうかを見てみましょう。

f:id:cross_hyou:20210923082758p:plain

year2019のp値は0.647、oldr:year2019のp値は0.511なので、年による違いがあるとは言えないようです。

2019年のときの回帰式は、

hosp = 237.2150 + 16.2060 + (4.1770 + 0.8714)*oldr

2007年のときの回帰式は、

hosp = 237.2150 + 4.1770*oldr

です。oldrの係数がyearを入れていないmodel1のときと比べるとだいぶ小さくなっていますね。

anova関数でmodel1とmodel2を比較してみましょう。

f:id:cross_hyou:20210923083403p:plain

p値がかなり小さい値ですので、model1とmodel2では有意な違いがあります。

つまり、model2でyear2019、oldr:year2019の個々の係数のp値は有意ではないですが、二つ合わせると有意ということですね。

anova関数だと簡単に計算できますが、これをマニュアルというか、なるべく原始的に計算するとどうなるでしょうか?

F = ( (SSRr - SSRur)/q ) / ( SSRur / (n - k - 1) )

というF値を算出して、F検定します。

SSRrはmdel1のほうの残差の2乗の合計(residual sum of squares)で、

SSRurはmodel2のほうの残差の合計です。

qはmodel2に追加した説明変数の数ですかたら、この場合はyear2019, oldr:year2019なので2になります。

nは観測数で、今回は47都道府県*2年なので94です。

kはmodel2の説明変数の数ですから、oldr, year2019, oldr:year2019の3になります。

まず、Fを計算します。

f:id:cross_hyou:20210923090045p:plain

F値が20.64868と算出されました。これは、anova関数の結果のFの値と同じですね。

pf関数でp値を計算します。

f:id:cross_hyou:20210923090536p:plain

4.162873e-08となりました。anova関数のp値と同じですね。(anova関数は4.163e-08)

散布図に回帰直線を入れてみましょう。

f:id:cross_hyou:20210923091313p:plain

f:id:cross_hyou:20210923091343p:plain

黒い直線がyearを説明変数に入れない回帰直線です。

赤い直線は2007年の、緑の直線が2019年の回帰直線です。

緑の直線の切片が赤い直線の切片よりも上にあります。2019年のほうが同じoldrだとしてもhospは大きいということですよね。

傾きも2007年よりも2019年のほうがきつくなっています。2019年のほうがoldrの影響がより大きくなっているということですね。

今回は以上です。

はじめから読むには

 

www.crosshyou.info

です。

都道府県別の通院者率のデータの分析4 - R言語で単回帰分析 - 65歳以上人口割合で通院者率を回帰分析

f:id:cross_hyou:20210922200840j:plain

Photo by Sana Ullah on Unsplash 

www.crosshyou.info

の続きです。

通院者率は、oldr: 65歳以上の人口割合がと関係があるかどうかを調べてみます。

やっぱり歳をとってくると病院のお世話になることが若いときよりも多くなると思うんですよね。

oldrがある年を確認します。

f:id:cross_hyou:20210922201653p:plain

2005年から2019年まで毎年、データがあるのですね。

一番古い年でoldrとhosp: 人口1000人当たりの通院者率のあるのは2007年、一番新しい年で両方にデータがあるのは2019年です。

この2007年と2019年のデータを使って、oldrとhospの散布図を描いてみます。

f:id:cross_hyou:20210922202403p:plain

f:id:cross_hyou:20210922202415p:plain

2007年と2019年で分布範囲は違っていますが、oldrとhospに正の相関関係があることは間違いないようですね。

lm()関数で回帰分析してみましょう。

2007年と2019年だけのデータフレームを作ります。

f:id:cross_hyou:20210922202906p:plain

worr: 従業者率とgymn: 人口100万人当たりの社会体育施設数は2007年と2019年はデータが無かったのですね。

yearというファクタ型の変数で2007と2019を区別できるようにしておきました。

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

f:id:cross_hyou:20210922203254p:plain

このモデル全体のp-valueは2.2e-16よりも小さいので有意な統計モデルです。

hosp = 156.4482 -+ 7.9904*oldr + u

という推定式です。uはerror termです。

つまり、oldrが1ポイント上昇すると、通院者率が約8人増えるということです。

残差プロットを描いてみます。

f:id:cross_hyou:20210922203643p:plain

f:id:cross_hyou:20210922203654p:plain

残差が適当に散らばってみえますので、heteroskedasticity(誤差項の不均一分散)ではないようです。

確認してみましょう。

 

にあるSpecial Case of the White Test for Heterosjedasticityの方法でやってみます。

f:id:cross_hyou:20210922205342p:plain

p-value: 0.3069ということなので、uhat2 = yhat + yhat^2 + error という回帰式は有意ではないです。つまり、heteroskedasticity(誤差項の不均一分散)ではないということですね。

今回は以上です。

次回は

 

www.crosshyou.info

です。

はじめから読むには

 

www.crosshyou.info

です。

 

都道府県別の通院者率のデータの分析3 - 沖縄県が一番、通院者率が低い。

f:id:cross_hyou:20210920193929j:plain

Photo by John Lee on Unsplash 

www.crosshyou.info

の続きです。

hosp: 人口1000人当たりの通院者率は何年に調査しているのかを確認します。

f:id:cross_hyou:20210920194336p:plain

1989年が一番古く、2019年が最新です。1995年と2016年は観測数が46とひとつ足りないです。

2019年では、どの都道府県別の通院者率が高い、低いでしょうか?

f:id:cross_hyou:20210920194548p:plain

f:id:cross_hyou:20210920194558p:plain

岩手県、秋田県、島根県が通院者率が高く、沖縄県、石川県、滋賀県が低いです。

1989年はどうでしょうか?

f:id:cross_hyou:20210920194812p:plain

f:id:cross_hyou:20210920194822p:plain

山形県、岩手県、広島県が高く、沖縄県、宮崎県、山梨県が低いです。

沖縄県は1989年も2019年も一番、通院者率が低いですね。

1989年と2019年の差を計算してみます。

f:id:cross_hyou:20210920195846p:plain

f:id:cross_hyou:20210920195858p:plain

青森県、島根県、山口県が大きく上昇、広島県、岡山県、愛知県があまり上昇していない(といっても100人以上増えていますが)ことがわかりました。

今回は以上です。

次回は

 

www.crosshyou.info

です。

はじめから読むには

 

www.crosshyou.info

です。

都道府県別の通院者率のデータの分析2 - 各変数の経年変化を箱ひげ図に原データを重ねて表現する。

f:id:cross_hyou:20210920134417j:plain

Photo by Julia Zolotova on Unsplash 

www.crosshyou.info

の続きです。

今回は各変数の経年変化をグラフにしてみてみます。

まずは、hosp: 人口千人当たりの通院者率のデータです。

まず、summary関数で平均値などを確認します。

f:id:cross_hyou:20210920134626p:plain

平均値は327.7人です。最小値は172.3人、最大値は461.7人です。

それでは経年変化をグラフにしてみてみましょう。

f:id:cross_hyou:20210920134839p:plain

f:id:cross_hyou:20210920134849p:plain

あきらかな上昇傾向ですね。

次は、oldr: 65歳以上の人口の割合(%)です。

サマリは以下のようになります。

f:id:cross_hyou:20210920135016p:plain

平均は26.11%、最小値は16.10%、最大値は37.20%です。

経年変化をグラフにします。

f:id:cross_hyou:20210920135143p:plain

f:id:cross_hyou:20210920135153p:plain

こちらもあきらかな上昇トレンドですね。日本の高齢化がはっきりわかります。

 

worr: 従業者の割合(%)はどうでしょうか?

f:id:cross_hyou:20210920135326p:plain

平均値は94.70%、最大値は88.10%、最大値は97.10%です。

経年変化はどうでしょうか?

f:id:cross_hyou:20210920135450p:plain

f:id:cross_hyou:20210920135500p:plain

4年間しかデータが無いですが、上昇傾向ではないですね。

 

gymn: 人口100万人当たりの社会体育施設数を見ましょう。

f:id:cross_hyou:20210920135620p:plain

平均値は432.1、最小値は25.9、最大値は2042.7です。

経年変化を見ます。

f:id:cross_hyou:20210920135759p:plain

f:id:cross_hyou:20210920135809p:plain

2000年以前はだんだんと増加していったのが2000年以降は横ばいですね。

そして、最大値を記録している年が1980年ごろというのが不思議です。2000以上あったのが急に減ったんでしょうか?

今回は以上です。

次回は、

 

www.crosshyou.info

です。

はじめから読むには

 

www.crosshyou.info

です。

 

都道府県別の通院者率のデータの分析1 - R言語にデータを読み込む。

f:id:cross_hyou:20210920081959j:plain

Photo by Ana Markovych on Unsplash 

今回は都道府県別の通院者率のデータを分析してみようと思います。

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

f:id:cross_hyou:20210920082143p:plain

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

f:id:cross_hyou:20210920082213p:plain

続いて、65歳以上人口割合、就業者率、社会体育施設数、通院者率の4項目をデータとして選択しました。

f:id:cross_hyou:20210920082309p:plain

ダウンロードしたCSVファイルはこのようなものです。10行目に変数名を作って挿入しています。

このCSVファイルのデータをR言語に読み込みます。

まず。tidyverseパッケージの読み込みをしておきます。

f:id:cross_hyou:20210920082624p:plain

read_csv関数で読み込みます。

f:id:cross_hyou:20210920082919p:plain

head関数でどんな感じで読み込まれたか見てみましょう。

f:id:cross_hyou:20210920083132p:plain

あ、今回は文字化けせずに上手く読み込んでいますね。

変数名を説明します。

ycod: year code - 調査年のコード

year: year - 調査年

pcod: prefecture code - 都道府県コード

pref: prefecture - 都道府県名

oldr: old ratio - 65歳以上の人の人口割合(%)

worr: work ratio - 就業者の人口割合(%)

gymn: gym number - 人口100万人当たりの社会体育施設数(数)

hosp: hospital - 人口1000人当たりの通院者率(人)

です。

ycodを100000を引いて1000000で割って4桁の西暦の数値に直しましょう。

yearとprefを文字列型からファクター型に変更しましょう。

f:id:cross_hyou:20210920084724p:plain

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

f:id:cross_hyou:20210920084940p:plain

あらら。。。

yearとprefが文字化けしてしまいました。。

仕方ないですね。yearとprefは削除してしまいます。

f:id:cross_hyou:20210920085149p:plain

pcodが1000は北海道で、47000は沖縄県なのですが、このままではわかりにくいので、あらかじめ用意してある、

f:id:cross_hyou:20210920090038p:plain

このCSVファイルを読み込んで、このデータと結合します。

f:id:cross_hyou:20210920090255p:plain

inner_join関数で結合します。dfのpcodとpref_codeのcodeが同じデータです。

f:id:cross_hyou:20210920090534p:plain

うまく結合できました。pcodはもう必要ないので削除して、変数の並び順を少し直します。

f:id:cross_hyou:20210920090806p:plain

できました。新しく追加された変数を説明します。

east: 東日本なら1、西日本なら0のダミー変数

big6: 東京都、千葉県、神奈川県、埼玉県、愛知県、大阪府なら1、その他は0のダミー変数

nose: 海が無い県は1、ある県は0のダミー変数

今回は以上です。

次回は

 

www.crosshyou.info

です。

 

OECD Meat Consumption Data Analysis 8 - Serial Correlation Robust Inference using R

f:id:cross_hyou:20210919201421j:plain

Photo by T o T on Unsplash 

www.crosshyou.info

This post is following of above post.

In this post, I will check if there is serial correlation in the previous regression model.

First, I make residual with resid() function.

f:id:cross_hyou:20210919201617p:plain

Then, let's make a graph for the residual.

f:id:cross_hyou:20210919201730p:plain

f:id:cross_hyou:20210919201740p:plain

I cannot tell if there is serial correlation by watching above graph.

Let's check it.

f:id:cross_hyou:20210919201912p:plain

st_residual and L(st_residual) are highly correlated. So, st_model has serial correlation.

So, let's make Serial Correlation Robust Inference with lmtest package and sandwich package.

f:id:cross_hyou:20210919202127p:plain

capi is still statistically significant, t-value is 9.2124.

Let's see usual inference

f:id:cross_hyou:20210919202315p:plain

capi's t-value is 10.859, it is greater than serial correlation robust inference.

That's it. Thank you!

To read the 1st post,

 

www.crosshyou.info

 

OECD Meat Consumption Data Analysis 7 - Time Series Regression using R dynlm() function.

f:id:cross_hyou:20210919153643j:plain

Photo by Ashutosh Saraswat on Unsplash 

www.crosshyou.info

This post is following of aabove post.
In this post, I will do some time-series regression with R.

First, I made JPN only dataframe.

f:id:cross_hyou:20210919154056p:plain

Let's see df_jpn.

f:id:cross_hyou:20210919154413p:plain

Then, I make ts object form df_jpn.

f:id:cross_hyou:20210919154507p:plain

f:id:cross_hyou:20210919154638p:plain

Let's see pokg and capi.

f:id:cross_hyou:20210919155529p:plain

f:id:cross_hyou:20210919155539p:plain

The both, pokg: POULTRY KG_CAP and capi: per capita GDP have uptrend.

I will make 4 time-series regression models.

1. Static Model

2. Finite Distributed Lag Model

3. Static Model with Trend

4. Finite Distrubuted Lag Model with Trend.

Before make regression models, I load dynlm package.

f:id:cross_hyou:20210919155958p:plain

All right, let's make those models.

f:id:cross_hyou:20210919160614p:plain

I use stargazer package to make results table.

f:id:cross_hyou:20210919160752p:plain

f:id:cross_hyou:20210919161001p:plain

We see capi has positive coefficients.

That's it. Thank you!

Next post is

 

www.crosshyou.info

 

To read the 1st post,

 

www.crosshyou.info