Photo by olena ivanova on Unsplash
の続きです。
前回の分析では、人口当たりのあんま・マッサージ師の数が大きく増えているところもあれば、大きく減少しているところもありました。
そこで、この人口当たりのあんま・マッサージ師の変化幅を被説明変数にして回帰分析をしてみます。
最初に確認として、あんま・マッサージ師の変化幅をヒストグラムにしてみます。
hist()関数を使います。
hist()関数でfreq = FALSE にすると、度数でなくて比率でヒストグラムを描きますので、density()関数とlines()関数で上のように赤い確率密度分布を追加できます。
-1500から1000と増加している県もあれば減っている県もあるのだとわかります。
説明変数は何がいいでしょうか？
まずは、このあんま・マッサージ師の変化幅と一番相関の強い変数を説明変数にしてみましょう。
相関係数を調べます。
はじめに相関係数を格納する変数として、correlationsという名前の変数を作っておきます。そして、for()関数でdf_7518の2列目から順番に相関係数を計算して、correlationsに保存しています。このとき変数名も一緒に保存しています。
harikyu_pop_chgが0.67で一番相関係数が高い、となりました。実際に確認してみます。
0.671755とfor()関数で計算した値と一致しました。
散布図を描いてみます。
右肩上がりの散布図ですね。人口100万人当たりのはり・きゅう師の数が増えているところは、あんま・マッサージ師も増えている、という関係ですね。
それでは、lm()関数で回帰分析をしてみます。
回帰分析での推計結果は、
あんま・マッサージ師の変化幅 = -204 + 0.383 * はり・きゅう師の変化幅 + 誤差項
になります。
つまり、はり・きゅう師の変化幅が100人だったら、あんま・マッサージ師の変化幅は38人ということですね。
summary()関数を使うともう少し詳しい結果が表示されます。
harikyu_pop_chgの標準誤差(Std. Error)は0.06292でt値は6.083で、p値(Pr(>|t|)は2.35e-07ですので、ほぼ0です。つまり、harikyu_pop_chgの係数 = 0 という帰無仮説は棄却できます。
broomパッケージのtidy()関数を使うと、回帰分析の結果をデータフレームにして出力してくれて、さらに信頼区間も表示してくれます。やってみます。
harikyu_pop_chgの係数(estimate)の95%の信頼区間は、0.256から0.509だとわかります。
回帰分析で気を付けなくてはいけないのは、誤差項の分散が均一かどうか(Homoskedasticity)です。誤差項をプロットしてみます。
なんか大丈夫そうな感じですね。
誤差項が不均一分散(Heteroskedasticity)になっていないかどうか、Breusch - Pagan 検定をしてみます。
lmtestというパッケージのbptest()関数で実行できます。
p値は0.66と0.05よりも大きな値なので、誤差項が均一分散しているという帰無仮説を棄却できないということです。なので、この回帰分析モデルは誤差項が均一分散だとみて大丈夫です。
最後に、実際のannmassage_pop_chgと回帰分析の推計式から予測された値を散布図にしてみます。
和歌山県は実際は200人ぐらい増えていますが、予測では1000人ぐらい減る予測ですね。
赤い線が 実際の値 = 予測の値 の線なので、この線に近いほど予測と実際が近いということです。
今回は以上です。
次回は
です。
初めから読むには、
です。
