crosshyou

主にクロス表(分割表)分析をしようかなと思いはじめましたが、あまりクロス表の分析はできず。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

です。

読書記録 - 「南極の氷に何が起きているか - 気候変動と氷床の科学」杉山 慎 著 中公新書

南極にある氷は、地球にある淡水の大部分だそうです。

そして、その氷が地球温暖化の影響で融けて少なくなってきているそうです。

ただ、南極の氷全体から見たら融けている量はわずからしいですが、そのわずかでも地球には大きな影響のようです。

 

都道府県別のあんま・マッサージ師、はり・きゅう師、柔道整復師数のデータの分析3 - 1975年のデータと2018年のデータの比較。

f:id:cross_hyou:20220123075050j:plain

Photo by Dawid Zawiła on Unsplash 

www.crosshyou.info

の続きです。

2018年のデータフレームと1975年のデータフレームを合体させましょう。

f:id:cross_hyou:20220123080310p:plain

まず、上のようにして、df_1975、df_2019のそれぞれのデータフレームから必要な変数だけを抜き出しました。今回はpop: 人口と、anmassage_pop: 人口100万人当たりのあんま・マッサージ師数、harikyu_pop: 人口100万人当たりのはり・きゅう師数、judo_pop: 人口100万人当たりの柔道整復師数を選択しました。そして、それぞれの変数名を1975年なのか2018年なのかを区別できるように名前を書き替えました。

この2つの~~~~_tempデータフレームをinner_join()関数でprefを鍵にして結合します。

f:id:cross_hyou:20220123080741p:plain

summary()関数でdf_7518を見てみましょう。

f:id:cross_hyou:20220123081009p:plain

ちょうどよく、pop, anmassage_pop, harikyu_pop, judo_popの1975年と2018年が上下に並んでくれました。人口以外はあきらかに2018年のほうが数値が増えていますね。

人口について、分布を比較してみましょう。

f:id:cross_hyou:20220123082228p:plain

f:id:cross_hyou:20220123082244p:plain

上のヒストグラムが1975年の人口のヒストグラムで、下のヒストグラムが2018年のヒストグラムです。あまり変化はないですね。

t.test()関数でt検定してみます。

f:id:cross_hyou:20220123082843p:plain

p値は0.005なので、両者の平均値は同じ、という帰無仮説を棄却します。

pop75の平均値は2,381,695人でpop18の平均値は2,690,277人ですので2018年のほうが人口の平均値は大きいのですね。

それぞれの変数の相関関係をみてみましょう。

f:id:cross_hyou:20220123083706p:plain

pop75とpop18の相関係数は0.984でかなり高いです。popが絡んでいない中では、anmassage_pop75とharikyu_pop75が0.795で一番高いです。同じ年の組み合わせを除くとannmassage_pop75とanmassage_pop18が0.412で一番高いです。

散布図のマトリックスも作成してみましょう。

f:id:cross_hyou:20220123084454p:plain

f:id:cross_hyou:20220123084507p:plain

変数が多すぎてよくわからないです。

df_7518のサマリーを見返すと、judo_popが一番変化率が大きいようです。

judo_pop75の平均値は83.6で、judo_pop18の平均値は504.9と6倍くらいに増えています。

どの都道府県が柔道整復師数の数を増やしているのか見てみましょう。

f:id:cross_hyou:20220123090044p:plain

f:id:cross_hyou:20220123085653p:plain

f:id:cross_hyou:20220123085711p:plain

大阪府が一番増えています。京都府、和歌山県と続いています。

あんま・マッサージ師の変化幅はどうでしょうか?

f:id:cross_hyou:20220123090435p:plain

f:id:cross_hyou:20220123090446p:plain

あんま・マッサージ師は減っているとkころもありますね。和歌山県、香川県などは減っています。神奈川県が一番多く増えていて、京都府、東京都と続きます。

はり・きゅう師はどうでしょうか?

f:id:cross_hyou:20220123091207p:plain

f:id:cross_hyou:20220123091217p:plain

はり・きゅう師も減少している都道府県があります。鳥取県と山口県です。

大阪府が一番増えて、2番目が京都府、3番目が神奈川県です。

今回は以上です。

次回は、

 

www.crosshyou.info

です。

初めから見るには、

 

www.crosshyou.info

です。

 

RのHistDataパッケージのBowley

f:id:cross_hyou:20220122190800j:plain

Photo by Greg Rosenke on Unsplash 

HistDataパッケージのBowleyデータをみてみます。

これは、一番初めに書かれた統計の教科書の一つ、Arthur Bowley(1901)のなかで使われているデータで、イギリスとアイルランドの輸出データの時系列データです。

Arthur Bowleyはこの時系列データをグラフにして、移動平均線を表示してるそうです。

データを呼び出してみます。

f:id:cross_hyou:20220122191226p:plain

str()関数でBowleyのデータ構造を調べたところ、45行、2列のデータフレームで、
Yearという名前の変数と、Valueという名前の変数があることがわかります。

Yearは西暦で、Valueが輸出の数量だか金額ですね。

早速、ヘルプに書かれているコードを実行してみましょう。

f:id:cross_hyou:20220122195347p:plain

f:id:cross_hyou:20220122195358p:plain

plot()関数で時系列のグラフを描きました。

続いてのコードはこれです。

f:id:cross_hyou:20220122195536p:plain

statsパッケージのfilter()関数で移動平均の値を計算します。runningという関数を作って、3年、5年、9年移動平均を作りました。

stats::filter(x, rep(1/width, width), sides = 2) というので、widthの値で移動平均の計算周期を調整するようです。

そして、この後、lines()関数で移動平均を描いたグラフに追加していきます。

f:id:cross_hyou:20220122195945p:plain

f:id:cross_hyou:20220122195959p:plain

移動平均線が追加されました。

この後、lowess()関数で平滑化した線を追加しています。

f:id:cross_hyou:20220122200120p:plain

f:id:cross_hyou:20220122200129p:plain

最後はggplot2()パッケージでグラフを描いています。

f:id:cross_hyou:20220122200318p:plain

f:id:cross_hyou:20220122200336p:plain

これでヘルプのコードは以上です。

以下が今回のコードです。

#
# HistDataの読み込み
library(HistData)
#
# Bowleyの呼び出し
data(Bowley)
#
# str()
str(Bowley)
#
# plot the data
with(Bowley, plot(Year, Value, type = "b", lwd = 2,
                  ylab = "Value of British and Irish Exports",
                  main = "Bowley's example of the method 
                  of smoothing curves"))
#
# find moving averages - use center alignment (requires width = 0DD)
require(gtools, warn.conflicts = FALSE)
#
# simpler version using stats::filter
running <- function(x, width = 5) {
  as.vector(stats::filter(x, rep(1/width, width), sides = 2))
}
mav3 <- running(Bowley$Value, width = 3)
mav5 <- running(Bowley$Value, width = 5)
mav9 <- running(Bowley$Value, width = 9)
#
lines(Bowley$Year, mav3, col = "blue", lty = 2)
lines(Bowley$Year, mav5, col = "green3", lty = 3)
lines(Bowley$Year, mav9, col = "brown", lty = 4)
#
# add lowess smooth
lines(lowess(Bowley), col = "red", lwd = 2)
#
if (require("ggplot2", quietly = TRUE)) {
  ggplot(aes(x = Year, y = Value), data = Bowley) +
    geom_point() +
    geom_smooth(method = "loess", formula = y ~x)
}

都道府県別のあんま・マッサージ師、はり・きゅう師、柔道整復師数のデータの分析2 - 人口当たりの数の散布図をplot()関数で描く。

f:id:cross_hyou:20220122075451j:plain

Photo by Federico Di Dio photography on Unsplash 

www.crosshyou.info

の続きです。

2018年だけのデータフレームを作成します。

f:id:cross_hyou:20220122075957p:plain

このデータフレームを使って、anmassage_pop, harikyu_pop, judo_popの相関関係をみてみます。これらは順番に人口100万人当たりのあんま・マッサージ師、はり・きゅう師、柔道整復師数です。

f:id:cross_hyou:20220122080425p:plain

あんま・マッサージ師とはり・きゅう師の相関係数は0.78で一番高いです。

2番目がはり・きゅう師と柔道整復師数の相関係数で0.69です。

3番目があんま・マッサージ師と柔道整復師の相関係数で0.50です。

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

f:id:cross_hyou:20220122080927p:plain

f:id:cross_hyou:20220122080937p:plain

こんどは、それぞれの散布図で都道府県名を表示してみます。まずは、anmassage_popとharikyu_popからやってみます。

f:id:cross_hyou:20220122082004p:plain

f:id:cross_hyou:20220122082018p:plain

plot()関数でtype="n"にして空白のグラフを作り、text()関数で都道府県名をプロットしています。lm()関数で回帰直線を作り、それをabline()関数でグラフに描いています。

この赤い直線の上にあるほど、あんま・マッサージ師と比較してはり・きゅう師が多いことになります。大阪府、京都府、和歌山県、奈良県、兵庫県と関西の府県がはり・きゅう師が多いことがわかります。

次は、あんま・マッサージ師と柔道整復師数の散布図を同じようにして作ります。

f:id:cross_hyou:20220122082754p:plain

f:id:cross_hyou:20220122082805p:plain

大阪府、京都府、和歌山県、富山県などがあんま・マッサージ師に比較して柔道整復師の数が多いですね。大阪、京都府、和歌山県ははり・きゅう師や柔道整復師の数に比べてあんま・マッサージ師の数が少ない、と表現したほうがいいかもです。

3つ目の散布図は、はり・きゅう師と柔道整復師です。

f:id:cross_hyou:20220122083629p:plain

f:id:cross_hyou:20220122083643p:plain

大阪府ははり・きゅう師、柔道整復師ともに全国で一番、人口当たりの数が多いのですね。

いままでは2018年のデータでした。

これと同じことを一番古い1975年のデータでやってみましょう。今から約50年前のデータです。

まずは相関係数と散布図マトリックスをみてみます。

f:id:cross_hyou:20220122084631p:plain

f:id:cross_hyou:20220122084600p:plain

あんま・マッサージ師とはり・きゅう師の相関係数は0.79(2018年は0.78)

あんま・マッサージ師と柔道整復師の相関係数は0.45(2018年は0.50)

はり・きゅう師と柔道整復師の相関係数は0.18(2018年は0.69)

となりました。はり・きゅう師と柔道整復師の相関関係が大きく変化していますね。

グラフにして「見える化」してみましょう。

f:id:cross_hyou:20220122090526p:plain

f:id:cross_hyou:20220122090539p:plain

legend()関数で左上に凡例を表示しています。

緑が2018年で赤が1975年です。1975年は香川県が一番人口当たりのはり・きゅう師が多かったのですね。大阪府はその他大勢の中に入っていてどこにあるかわからないですね。全般的にはり・きゅう師、柔道整復師ともに増加しています。

回帰直線の傾きは2018年のほうが急勾配になっています。

今回は以上です。

次回は

 

www.crosshyou.info

です。

初めから読むには、

 

www.crosshyou.info

です。

都道府県別のあんま・マッサージ師、はり・きゅう師、柔道整復師数のデータの分析1 - 人口当たりでは、大阪がはり・きゅう師、柔道整復師数が一番多い。

f:id:cross_hyou:20220115192214j:plain

Photo by Alex Machado on Unsplash

今回は都道府県別のあんま・マッサージ師、はり・きゅう師、柔道整復師数のデータを分析してみようと思います。

まず、政府統計の総合窓口、e-stat.go.jpからデータをダウンロードします。

f:id:cross_hyou:20220115192359p:plain

47の都道府県を選択し、

f:id:cross_hyou:20220115192421p:plain

総人口、あんま・マッサージ師数、はり・きゅう師数、柔道整復師数を選択しました。

f:id:cross_hyou:20220115192505p:plain

このようなCSVファイルを取得できたので、これをRで分析します。

f:id:cross_hyou:20220115192809p:plain

前準備としてtidyverseパッケージの読み込みをしておきます。

read_csv()関数でCSVファイルを読み込みます。

f:id:cross_hyou:20220116082344p:plain

head()関数で頭の数行を表示してみます。

f:id:cross_hyou:20220116082555p:plain

yearを100000を引いてから1000000で割って普通の4桁の西暦に戻すのとNAのある行を削除しましょう。

f:id:cross_hyou:20220116083015p:plain

summary()関数でdfを見てみます。

f:id:cross_hyou:20220116083139p:plain

yearが4桁の西暦になっていることがわかります。一番古いデータは、1975年で、一番新しいデータは2018年です。

year2は195年度などの文字列です。

codeは北海道が1000, 沖縄が47000です。

prefは都道府県名の文字列です。

popは総人口、一番少ない都道府県・年は560,000人、一番多い都道府県・年は13,822,000人です。

anmassageはあんま・マッサージ師の数です。一番少ないのは61人、多いのは27,399人です。

harikyuはり・きゅう師の数です。一番少ないのは138人、多いのは43,909人です。

judoは柔道整復師数です。一番少ないのは3人、多いのは11,105人です。

はり・きゅう師 > あんま・マッサージ師 > 柔道整復師数 という感じの人数ですね。

年ごとの合計のデータフレームを作って、それぞれの経年変化をみてみます。

f:id:cross_hyou:20220116084405p:plain

group_by()関数やsummarize()関数を使って、年ごとの合計値のデータフレームを作成しました。これをplot()関数とline()関数でグラフにて、legend()関数で凡例を表示します。

f:id:cross_hyou:20220116085815p:plain

f:id:cross_hyou:20220116085826p:plain

はり・きゅう師の伸びが加速している感じですね。

今度は、都道府県別の状況を見てみましょう。

人口百万人当たりの数を計算して、それをくらべましょう。

f:id:cross_hyou:20220116181301p:plain

こんな感じでmutate()関数で人口当たりの数の変数を作りました。

一番新しい2018年のデータで比べてみます。

f:id:cross_hyou:20220116181918p:plain

f:id:cross_hyou:20220116181929p:plain

東京都がダントツに多いですね。人口当たりの数のであんま・マッサージ師は2000人近くいます。岩手県が一番少なく500人以下ですね。

はり・きゅう師はどうでしょうか?

f:id:cross_hyou:20220116182347p:plain

f:id:cross_hyou:20220116182357p:plain

はり・きゅう師は大阪が一番ですね。東京は2番です。

柔道整復師数はどうでしょうか?

f:id:cross_hyou:20220116182751p:plain

f:id:cross_hyou:20220116182801p:plain

大阪が一番多いですね。

今回は以上です。

次回は

 

www.crosshyou.info

です。

 

















 

RのHistDataパッケージのArmada

f:id:cross_hyou:20220115075815j:plain

Photo by Peter Pryharski on Unsplash 
RのHistDataパッケージのArmadaデータセットはスペインの無敵艦隊のデータです。

f:id:cross_hyou:20220115080128p:plain

str()関数でデータ構造を見ると、10行 x 11列のデータフレームです。

それぞれの変数が何かというと、

Armada: designation of the fleetとあります。艦隊の行先でしょうか?

ships: 船の数

tons: total tonsとあります。艦隊のトン数ですかね。。

soldiers: 兵士の数

sailors: 船員の数

men:  兵士の数と船員の数の合計

artillery: 大砲の数

balls: 玉の数です。砲弾の数ですかね?

gunpowder: 火薬の数

lead: リードの数、リードってなんでしょうかね?海軍というか船関係のものだと思いますが、わからないです。

roap: ロープの数、ロープは紐のロープでしょうね。。

ヘルプページにあるコードを動かしてみます。

まずは、単純に全データの表示です。

  

f:id:cross_hyou:20220115081640p:plain

ArmadaはPortugalやNapolesなど地名になっているので艦隊の行先なのですかね?

その他の変数はすべて数値データです。

続いてのヘルプのコードを実行します。

f:id:cross_hyou:20220115081921p:plain

armada <- Armada[ , -c(1, 6)] で1番目の変数、Armadaと6番目の変数、menを削除して新しいデータフレーム、armadaを作っています。men = soldiers + sailors なのでいらないということですね。次に、prcomp()関数でPCA分析ですね。Principal Component Analysys, 日本語だと主成分分析です。summary()関数で結果を表示しています。元のデータフレームは9つの変数がありますが、それをPC1, PC2 という合成変数にすると、この2つだけで87%ほど全部のデータのバラツキを説明できる、ということですね。

さらにヘルプのコードを実行します。

f:id:cross_hyou:20220115082704p:plain

f:id:cross_hyou:20220115082716p:plain

このコードは、PCAのサマリーの表の proportion of variance をグラフにしたものですね。PC1だけで全体の変動の7割以上を表現できるということですね。

続いてのヘルプのコードを実行します。

f:id:cross_hyou:20220115083008p:plain

f:id:cross_hyou:20220115083019p:plain

biplot()関数で、PC1とPC2のグラフにしています。

以上でヘルプのコードは終了です。

下に記載します。

library(HistData)
data("Armada")
str(Armada)
Armada
armada <- Armada[ , -c(1, 6)]
armada.pca <- prcomp(armada, scale = TRUE)
summary(armada.pca)
plot(armada.pca, type = "lines", pch = 16, cex = 2)
biplot(armada.pca)