Rで何かをしたり、読書をするブログ

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

都道府県別のあんま・マッサージ師、はり・きゅう師、柔道整復師数のデータの分析4 - Rのlm()関数で単回帰分析をする。broomパッケージやlmtestパッケージも使用。

f:id:cross_hyou:20220129083625j:plain

Photo by olena ivanova on Unsplash 

www.crosshyou.info

の続きです。

前回の分析では、人口当たりのあんま・マッサージ師の数が大きく増えているところもあれば、大きく減少しているところもありました。

そこで、この人口当たりのあんま・マッサージ師の変化幅を被説明変数にして回帰分析をしてみます。

最初に確認として、あんま・マッサージ師の変化幅をヒストグラムにしてみます。

hist()関数を使います。

f:id:cross_hyou:20220129084757p:plain

f:id:cross_hyou:20220129084810p:plain

hist()関数でfreq = FALSE にすると、度数でなくて比率でヒストグラムを描きますので、density()関数とlines()関数で上のように赤い確率密度分布を追加できます。

-1500から1000と増加している県もあれば減っている県もあるのだとわかります。

説明変数は何がいいでしょうか?

まずは、このあんま・マッサージ師の変化幅と一番相関の強い変数を説明変数にしてみましょう。

相関係数を調べます。

f:id:cross_hyou:20220129090103p:plain

はじめに相関係数を格納する変数として、correlationsという名前の変数を作っておきます。そして、for()関数でdf_7518の2列目から順番に相関係数を計算して、correlationsに保存しています。このとき変数名も一緒に保存しています。

harikyu_pop_chgが0.67で一番相関係数が高い、となりました。実際に確認してみます。

f:id:cross_hyou:20220129090536p:plain

0.671755とfor()関数で計算した値と一致しました。

散布図を描いてみます。

f:id:cross_hyou:20220129090838p:plain

f:id:cross_hyou:20220129090849p:plain

右肩上がりの散布図ですね。人口100万人当たりのはり・きゅう師の数が増えているところは、あんま・マッサージ師も増えている、という関係ですね。

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

f:id:cross_hyou:20220129091230p:plain

回帰分析での推計結果は、

あんま・マッサージ師の変化幅 = -204 + 0.383 * はり・きゅう師の変化幅 + 誤差項

になります。

つまり、はり・きゅう師の変化幅が100人だったら、あんま・マッサージ師の変化幅は38人ということですね。

summary()関数を使うともう少し詳しい結果が表示されます。

f:id:cross_hyou:20220129091627p:plain

harikyu_pop_chgの標準誤差(Std. Error)は0.06292でt値は6.083で、p値(Pr(>|t|)は2.35e-07ですので、ほぼ0です。つまり、harikyu_pop_chgの係数 = 0 という帰無仮説は棄却できます。

broomパッケージのtidy()関数を使うと、回帰分析の結果をデータフレームにして出力してくれて、さらに信頼区間も表示してくれます。やってみます。

f:id:cross_hyou:20220129092222p:plain

harikyu_pop_chgの係数(estimate)の95%の信頼区間は、0.256から0.509だとわかります。

回帰分析で気を付けなくてはいけないのは、誤差項の分散が均一かどうか(Homoskedasticity)です。誤差項をプロットしてみます。

f:id:cross_hyou:20220129092610p:plain

f:id:cross_hyou:20220129092623p:plain

なんか大丈夫そうな感じですね。

誤差項が不均一分散(Heteroskedasticity)になっていないかどうか、Breusch - Pagan 検定をしてみます。

lmtestというパッケージのbptest()関数で実行できます。

f:id:cross_hyou:20220129093112p:plain

p値は0.66と0.05よりも大きな値なので、誤差項が均一分散しているという帰無仮説を棄却できないということです。なので、この回帰分析モデルは誤差項が均一分散だとみて大丈夫です。

最後に、実際のannmassage_pop_chgと回帰分析の推計式から予測された値を散布図にしてみます。

f:id:cross_hyou:20220129094747p:plain

f:id:cross_hyou:20220129094817p:plain

和歌山県は実際は200人ぐらい増えていますが、予測では1000人ぐらい減る予測ですね。

赤い線が 実際の値 = 予測の値 の線なので、この線に近いほど予測と実際が近いということです。

今回は以上です。

次回は

 

www.crosshyou.info

です。

初めから読むには、

 

www.crosshyou.info

です。