crosshyou

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

国税局別の民間給与実態調査のデータの分析6 - R言語で1~4人の事業所の人数と5,000人以上の事業所の人数を計算する。

 

www.crosshyou.info

 の続きです。

今回は、事業所の規模が1~4人の人数と5,000人以上の人数に注目してみます。

まず、spread関数で、1~4人の人数の列と5,000人以上の人数の列を持つデータフレームを作ります。

f:id:cross_hyou:20201123145921p:plain

変数名を変更します。rename関数を使います。

f:id:cross_hyou:20201123150227p:plain

smallが1~4人の事業所の人数で、largeが5,000人以上の事業所の人数です。

smallとlargeの比率も作ります。mutate関数です。

f:id:cross_hyou:20201123150545p:plain

sl_ratioが比率です。最小値は0.3131, 最大値は3.2966, 平均値は1.5602, 中央値は1.3015です。

geom_histogram関数でヒストグラムを描いてみます。

f:id:cross_hyou:20201123150907p:plain

f:id:cross_hyou:20201123150925p:plain

bin = 30だと細かすぎのようですね。bin = 10でやってみます。

f:id:cross_hyou:20201123151205p:plain

f:id:cross_hyou:20201123151222p:plain

1より小さいところと3より小さいところに頂点がある、2つの頂点の分布のようです。2015年と2016年で違うかyearendとaverageで違うか、facet_grid関数で処理してみます。

f:id:cross_hyou:20201123151820p:plain

f:id:cross_hyou:20201123151841p:plain

どうなんでしょうか。。。yearやitemによる分布の違いは無さそうですね。

lm関数でsl_ratioをresponse variable, yearとitemをexplanatory variablesにした回帰分析モデルを作ってみます。

f:id:cross_hyou:20201123152300p:plain

p-valueが0.9581ですから有意なモデルではないです。つまり、sl_ratioはyearやitemとは関連は無い、ということですね。

今回は以上です。

 

国税局別の民間給与実態調査のデータの分析5 - R言語で2x2のクロス表を分析する。chisq.test関数、fisher.test関数、assocstats関数、oddsratio関数。

 

www.crosshyou.info

 の続きです。

今回はR言語で2x2のクロス表を分析してみようと思います。

 

カテゴリカルデータ解析 (Rで学ぶデータサイエンス 1)
 

 こちらの本を参考にしました。

まず、table関数で2x2のクロス表を作ります。yearendとaverage, zouとgenでのクロス表を作ります。

f:id:cross_hyou:20201123091021p:plain

yearend(年末時点の人数)は増加が94個、減少が23個

average(年間平均の人数)は増加が90個、減少が27個だとわかります。

prop.table関数で比率であらわしてみます。

f:id:cross_hyou:20201123091413p:plain

yearendはzouが80%, genが20%, averageはzouが77%, genが23%とわかります。

このように、yearendとaverageではzouとgenの割合が有意に違うのか、カイ二乗検定で調べてみます。chisq.test関数を使います。

f:id:cross_hyou:20201123091805p:plain

p-value = 0.6323と0.05よりも大きな値です。つまり、yearendとaverageではzou, genの割合に有意な違いは無い、ということですね。

フィッシャーの直接確率法での検定もしてみます。fisher.test関数です。

f:id:cross_hyou:20201123092543p:plain

p-value = 0.6326 と0.05よりも大きく、カイ二乗検定と同じく有意な違いはありません。

vcdパッケージの中にassociate関数でファイ係数がわかります。これは2x2のクロス表の相関係数です。

f:id:cross_hyou:20201123093136p:plain

Phi-Coefficient : 0.042とありますので、ほとんど関連は無い、ということです。

Phi-Coefficientは2x2のクロス表では、相関係数と同じですから、相関係数を計算してみましょう。

f:id:cross_hyou:20201123094246p:plain

yerendを0, averageを1, zouを0, genを1に見立ててyearend, averageのベクトルとzou, genのベクトルを作り、cor関数で相関係数を計算しています。0.04170288とファイ係数と同じ値になりました。

oddsratio関数でオッズレシオを計算してみます。

f:id:cross_hyou:20201123094806p:plain

オッズレシオは1.226087になりました。

オッズレシオは2つのオッズの比ですから、手計算で求めてみましょう。

まず、yearendのzouのオッズです。これは、yearendのzou / yearendのgen で計算します。

f:id:cross_hyou:20201123095006p:plain

次にaverageのzouのオッズを計算します。

f:id:cross_hyou:20201123095106p:plain

この2つのオッズの比を計算すると、オッズレシオになります。

f:id:cross_hyou:20201123100035p:plain

oddsratio関数での結果と一致しました。

オッズレシオの95%信頼区間を求めましょう。

f:id:cross_hyou:20201123100705p:plain

95%の信頼区間が1.0を含んでいますので、やっぱりyearendとaverageには有意な違いは無いということです。

今回は以上です。

国税局別の民間給与実態調査のデータの分析4 - R言語で2015年と2016年の伸び率を計算する。

 

www.crosshyou.info

 の続きです。

今回は、2016年の人数が2015年と比較して、どれだけ伸びたかを計算したいと思います。

まず。df_longのデータフレーム、

f:id:cross_hyou:20201122164609p:plain

これを、

f:id:cross_hyou:20201122164632p:plain

こういうデータフレームに変換します。

spread関数を使います。

f:id:cross_hyou:20201122164851p:plain

変数名が2015年、2016年とはじめの文字が数字だとよくないので、修正します。

rename関数を使います。

f:id:cross_hyou:20201122165033p:plain

次はmutate関数を使って2015年から2016年の伸び率を作ります。

f:id:cross_hyou:20201122165230p:plain

officeが文字列型なので、as.factor関数でファクター型に変換しておきましょう。

f:id:cross_hyou:20201122165344p:plain

summary関数で、df_yearのサマリーを見てみます。

f:id:cross_hyou:20201122165503p:plain

chg_pctが伸び率です。最小値は-18.08%、最大値は56.32%、平均値は2.41%、中央値は1.76%です。前年から比べて、56%も増えるってすごいですね。

zougenのところ見ると、zou,増加したのが184個で、gen, 減少したのが50個あるとわかります。

それではchg_pctで並び替えてみます。arrange関数を使います。

f:id:cross_hyou:20201122170235p:plain

一番人数の伸び率が大きかったのは、5~9人の札幌国税局の年間平均人数です。

19万6460人が30万7097人に増えました。56.3%の伸び率です。

二番目が5,000人以上の熊本国税局の年末人数です。6万4663人が8万5121人にふえました。31.6%の伸び率です。

三番目が1~4人の仙台国税局の年間平均です。21万1152人が25万9064人に増えました。22.7%の伸びです。

 

減少したところはどこでしょうか?

f:id:cross_hyou:20201122170912p:plain

一番減少したのは、5~9人の熊本国税局の年末人数でした。24万5527人から20万1140人に、マイナス18.1%の減少です。

二番目に減少したのは、1~4人の福岡国税局の年間平均人数でした。22万5577人が18万5027人に、マイナス18.0%の減少です。

三番目に減少したのは、5~9人の熊本国税局の年間平均人数でした。24万4004人から20万2389人に、まいなす17.1%の減少です。

熊本国税局は増加でも、減少でも顔を出していますね。

itemはyearend, kiboは合計だけに絞ってみてみましょう。

f:id:cross_hyou:20201122171825p:plain

札幌国税局が189万2119人から196万8759人に、4.05%の上昇でいちばんでした。東京国税局、大阪国税局、全国合計と続きます。全ての国税局で増加していますね。関東信越国税局が一番伸びが小さいです。

今回は以上です。

 

国税局別の民間給与実態調査のデータの分析3 - R言語のgeom_bar関数で棒グラフを描く。function関数、lapply関数、grid.arrange関数で効率よく。

 

www.crosshyou.info

 の続きです。

今回はR言語のgeom_bar関数で棒グラフを描いてみます。

例えば、2015年、yearend, 1~4人、全国除くだと、

f:id:cross_hyou:20201122090635p:plain

f:id:cross_hyou:20201122090657p:plain

となります。tokyo, osaka, kantoshinetsu, nagoya, fukuokaという順番です。

facet_grid関数を使って、yearとitemの 2x2の棒グラフを描きます。

f:id:cross_hyou:20201122091436p:plain

f:id:cross_hyou:20201122091456p:plain

これを他の規模でも同じようにします。

まずは、function関数で上のコードを再現する関数を作ります。

f:id:cross_hyou:20201122093245p:plain

この関数をsapply関数を5~9人で試してみましょう。

f:id:cross_hyou:20201122093746p:plain

f:id:cross_hyou:20201122093802p:plain

うまくできました。

lapply関数で、全てのkiboについて、一つのコマンドで実行します。

まず、gridExtraパッケージを読み込みます。

f:id:cross_hyou:20201122095422p:plain

そして、lapply関数で各kiboのグラフオブジェクトを作成します。

f:id:cross_hyou:20201122095613p:plain

このkibo_graphがうまくいったか確認してみます。

f:id:cross_hyou:20201122095934p:plain

f:id:cross_hyou:20201122095952p:plain

うまくいきました。

それでは、gridExpandパッケージのgrid.arrange関数で9つのグラフをいっぺんに表示しましょう。

f:id:cross_hyou:20201122100448p:plain

f:id:cross_hyou:20201122100737p:plain

このようになります。事業所の規模が大きいところは、tokyoの人数が突出していますね。

今回は以上です。

 

国税局別の民間給与実態調査のデータの分析2 - R言語でデータフレームの形態を変更(横長型から縦長型に) - gather関数とspread関数

 

www.crosshyou.info

の続きです。

前回作成したデータフレーム, dfのはじめの数行をhead関数でみてみましょう。

f:id:cross_hyou:20201121180708p:plain

このように人数という値が何列にも配置されています。

横長型といえばいいのかな。。。

これを縦長型に変更してみたいと思います。

図で説明すると、

f:id:cross_hyou:20201121180931p:plain

となっているデータフレームに、下の図のように、

officeとninzuuいう新しい変数を作って、

officeにはsapporo, sendai, kantoshinetsuなどの変数名を格納し、

ninzuuには2083522などの値を格納します。

f:id:cross_hyou:20201121181308p:plain

tidyverseパッケージを読み込むとgather関数が利用できます。

この関数を使います。

まず、library関数でtidyverseを読み込みます。

f:id:cross_hyou:20201121182430p:plain

gather関数を使います。

gather(データフレーム, key = 元の変数名が入る新しい列名, value = 元のデータフレーム値が入る新しい列名, -valueに入れない変数名)

こんな構文です。

f:id:cross_hyou:20201121182905p:plain

gather関数は横長のデータフレームを縦長にします。

その反対がspread関数です。

spread(データフレーム, key = 新しく列名にしたい要素が入っている変数, value = keyに格納する要素が入っている変数)

という構文です。使ってみます。

f:id:cross_hyou:20201121183506p:plain

このgather関数とspread関数、正直言ってよくわかっていません。毎回、試行錯誤です。

今回は以上です。





 

国税局別の民間給与実態調査のデータの分析1 - R言語でデータを読み込む。

今回は国税局別の民間給与実態調査のデータを分析してみます。

いつものように、政府統計の総合窓口(www.e-stat.go.jp)からデータを取りました。

f:id:cross_hyou:20201121090804p:plain

f:id:cross_hyou:20201121090933p:plain

データベースのほうをクリックしました。

f:id:cross_hyou:20201121091051p:plain

年次をクリックします。

f:id:cross_hyou:20201121091200p:plain

2016年をクリックします。

f:id:cross_hyou:20201121091523p:plain

たくさん種類がありましたが、今回は【参考】国税局別表 第1表 を選択しました。DBのところをクリックします。

f:id:cross_hyou:20201121091752p:plain

こういうデータでした。人数のデータのようです。ダウンロードをクリックします。

f:id:cross_hyou:20201121092003p:plain

このようなポップアップが出てきました。ダウンロードをクリックします。

f:id:cross_hyou:20201121092114p:plain

さらにダウンロードをクリックします。

f:id:cross_hyou:20201121092338p:plain

こういうデータファイルでした。

これをcsvファイルとして保存してR言語に読み込みます。

f:id:cross_hyou:20201121092823p:plain

データが10行目から始まっているので、skip = 9 として始めの9行は無視しています。

データが無いところは、***か-なので、na.strings = c("***", "-")として***と-はNAにするようにしています。ファイルのエンコードがUTF-8なので、encoding = "UTF-8"としています。

str関数でファイルが読み込まれているかどうか見てみます。

f:id:cross_hyou:20201121093306p:plain

うまく読み込まれたようです。変数名を修正しましょう。

f:id:cross_hyou:20201121093856p:plain

もう一度、str関数で確認します。

f:id:cross_hyou:20201121094105p:plain

変数名が変わりました。

year, item, kibo, xは数値データではないので、ファクター型に変換します。

f:id:cross_hyou:20201121094900p:plain

as.factor関数でファクター型に変換しています。

また、str関数でファクター型に変換されているかみてみます。

f:id:cross_hyou:20201121095046p:plain

x は全部NAなので、必要ないです。削除します。

f:id:cross_hyou:20201121095332p:plain

str関数で見てみます。

f:id:cross_hyou:20201121095500p:plain

xがなくなっています。

yearは2015年と2016年の2つですね。

itemは2つですが、長い文字列なので、全部でてないですね。

levels関数で確認してみます。

f:id:cross_hyou:20201121095728p:plain

12月末の人数と年間平均の人数です。yearendとaverageに変更しましょう。

f:id:cross_hyou:20201121100050p:plain

table関数で各要素が何個あるか数えました。両方とも18個ありました。

kiboのlevelも確認します。

f:id:cross_hyou:20201121100403p:plain

これは事業所の人数の規模です。ordered関数でlevelの順番を変更しましょう。

f:id:cross_hyou:20201121101345p:plain

summary関数で各データのサマリーを見てみます。

f:id:cross_hyou:20201121101702p:plain

f:id:cross_hyou:20201121101642p:plain

これで分析のためのデータフレームが整いました。

今回は以上です。

 

都道府県別別の「医療費の動向」調査のデータ分析6 - R言語のglm関数でロジスティクス回帰分析をする。

 

www.crosshyou.info

 の続きです。

今回はR言語のglm関数でロジスティクス回帰分析をします。

response variableはeastweat: 東日本と西日本にしてみます。

explanatory variablesはyear: 年度、medical_total: 医科_計、dental: 歯科、pharma: 調剤、visit: 訪問看護療養にしてみます。

まずは、それぞれのexplanatory variableとeastjapanの関係を箱ひげ図であらわしてみます。

まずは、medical_totalから。

f:id:cross_hyou:20201120191824p:plain

f:id:cross_hyou:20201120191838p:plain

ggplot2パッケージのgeom_boxplot関数で箱ひげ図を描きました。

東日本と西日本では明確な違いは無いようですね。

次はdental: 歯科です。

f:id:cross_hyou:20201120192153p:plain

f:id:cross_hyou:20201120192208p:plain

geom_boxplot(aes(colour = eastwest))として東日本と西日本で色を変えてみました。

次はpahrma: 調剤です。

f:id:cross_hyou:20201120204459p:plain

f:id:cross_hyou:20201120204514p:plain

geom_point関数で個々のデータもプロットしました。

visit: 訪問看護療養はどうでしょうか?

f:id:cross_hyou:20201120205048p:plain

f:id:cross_hyou:20201120205103p:plain

geom_boxplot(aes(fill = eastjapan))として色の分けしたのと、

geom_point(aes(colour = year))として、個々のデータを色分けしてみました。

こうして箱ひげ図にしてみましたが、東日本と西日本ではっきりした違いはなさそうですね。

続いて、explanatory variablesどうしの相関係数マトリックスを見てみましょう。

GGallyパッケージのggpairs関数を使います。

f:id:cross_hyou:20201120205558p:plain

f:id:cross_hyou:20201120205526p:plain

一番相関係数の大きいのはmedical_totalとpharmaです。0.691です。explanatory variablesどうしの相関は高くないので、多重共線性はないですね。

それでは、glm関数でロジスティクス回帰分析のモデルを作ります。

f:id:cross_hyou:20201120210045p:plain

medical_total: 医科_計、dental: 歯科、visit: 訪問看護療養が有意な変数のようですね。

predict関数でこのモデルから東日本と西日本を予測してみます。

f:id:cross_hyou:20201120210823p:plain

table関数で予測と実際の2x2のクロス表を作ります。

f:id:cross_hyou:20201120211133p:plain

1が西日本、2が東日本です。

モデルが西日本と予測したうち、77は正解、31は不正解

モデルが東日本と予測したうち、12は正解、20は不正解です。

chisq.test関数でこの予測結果が偶然でもなるのか、統計的に有意なのか調べます。

f:id:cross_hyou:20201120211423p:plain

p-value = 0.001742と0.05より小さいので有意なモデルです。

オッズレシオを計算してみましょう。

f:id:cross_hyou:20201120212920p:plain

オッズレシオは分子が正解の数の掛け算、分母が不正解の数の掛け算です。

オッズレシオは3.82です。

今回は以上です。