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

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

商業統計調査データの分析4 - R言語で法人と個人の事業所の割合を比較する。

 

www.crosshyou.info

 の続きです。

今回は、R言語で法人の事業所の割合と個人の事業所の割合を出して、経年変化をみてみたいと思います。

まずは、もとになるデータをhead関数とsummary関数で確認しましょう。

f:id:cross_hyou:20190406145917j:plain

f:id:cross_hyou:20190406145929j:plain

このようなデータです。Corp_Totalが法人と個人を合わせた事業所の数、Corp_Corpが法人の事業所の数、Corp_Indiが個人の事業所の数でした。

Corp_Corp / Corp_Total で法人の割合、Corp_Indi / Corp_Totalで個人の割合になります。早速作成しましょう。

f:id:cross_hyou:20190406150706j:plain

まずは法人の割合を作成しました。平均値と中央値が同じ値ですね。0.4082です。最大値は0.8766、最小値は0.0244です。

個人の割合はこうなりました。

f:id:cross_hyou:20190406150834j:plain

個人の割合も平均値と中央値が同じになりました。

こうしてみると、だいたい6割が個人の事業所、4割が法人の事業所だとわかりました。

hist関数でヒストグラムを作成してみます。

f:id:cross_hyou:20190406151213j:plain

f:id:cross_hyou:20190406151226j:plain

前回の一人当り売上高、売場面積当り売上高と比べると左右対称の度合いが高いですね。

年毎の比率を見てみましょう。1 - 法人の割合 = 個人の割合ですから、ここからは法人の割合だけを見ていきます。tapply関数を使います。

f:id:cross_hyou:20190406151614j:plain

1988年は0.32ぐらいでしたが、年が経つたびに上昇し、2012年には0.499と半分近くにまで高まっています。

平均値だけでなくて、年ごとの箱ひげ図を作成します。Yearをfactor関数でファクタにしてX軸にしてplot関数です。

f:id:cross_hyou:20190406152048j:plain

f:id:cross_hyou:20190406152100j:plain

法人の割合 = a + b * 年
という線形モデルを作成して年が法人の割合に有意に影響しているかを調べましょう。

lm関数を使います。

f:id:cross_hyou:20190406152421j:plain

p-value: < 2.2e-16とありますので、このモデルは有意です。Yearの行のPr(>|t|)を見ると、<2e-16とありますので、Yearは有意な影響を及ぼしています。

plot関数で散布図を描いて、abline関数でモデルの回帰直線を重ねてみます。

f:id:cross_hyou:20190406152933j:plain

f:id:cross_hyou:20190406152945j:plain

こんな感じです。モデルの残差をグラフにしてみましょう。

f:id:cross_hyou:20190406153206j:plain

f:id:cross_hyou:20190406153218j:plain

このようになりました。

ところで、このYearだけで法人の割合を回帰するmodel1は、Adjusted R-squaredが0.1301とかなり低いです。他の変数も加えてもう少しAdjusted R-squaredの数値をあげてみたいです。法人の割合が多いのはどんな場合でしょうか?一事業所当りの売上高が多い業種は法人の割合が多いような気がします。やってみましょう。まず、一事業所当りの売上高の変数を作成します。

f:id:cross_hyou:20190406153828j:plain

一事業所当りの売上高の平均値は役1億2100万円です。中央値は約5600万円です。

この変数を加えて回帰分析をしてみましょう。

f:id:cross_hyou:20190406154239j:plain

一番下の行のp-valueは2.2e-16よりも小さいの有意です。Yearの係数もRev_Corの係数もp値は2e-16以下なので有意です。Adjusted R-Squaredは0.2554ですので、model1よりも高くなりました。

回帰残差をプロットしてみます。

f:id:cross_hyou:20190406154917j:plain

こうなりました。

anova関数を使って、model1とmodel2が有意に異なるかチェックします。

f:id:cross_hyou:20190406155029j:plain

p値が2.2e-16以下なので有意に異なります。model1に一事業所当り売上高を加えたことは意味があったということです。

一事業所当りの従業者数も加えてみましょう。一事業者当りの従業者数が大きいということは大きな事業所ということを示唆しますし、大きな事業所ということは法人の運営する事業所だと思われます。

f:id:cross_hyou:20190406155527j:plain

平均すると一事業所当り約6人いますね。

それではこのStf_Corを加えたmodel3を作成してみましょう。

f:id:cross_hyou:20190406155926j:plain

モデル全体はp-value* < 2.2e-16とありますので有意ですが、Rev_Cor, Stf_Corの係数が有意ではないですね。Rev_Cor(一事業所当り売上高)を除いて、Stf_Cor(一事業所当り従業者数)だけを残したmodel4を試してみます。

f:id:cross_hyou:20190406160406j:plain

model4は、Adjusted R-squaredが0.2605と一番いいです。

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

f:id:cross_hyou:20190406160709j:plain

f:id:cross_hyou:20190406160723j:plain

こうなりました。

回帰分析はここまでにして、法人の割合の高い業種と個人の割合の高い業種を調べてみましょう。tapply関数とsort関数とrev関数とhead関数です。

f:id:cross_hyou:20190406161054j:plain

法人の割合の高い業種は上の業種でした。各種商品小売業が約8割です。以下、燃料小売業、その他の各種商品小売業(従業者が常時50人未満のもの)、楽器小売業、農業用機械器具小売業、その他の機械器具小売業と続きます。

個人の割合の高い業種は

f:id:cross_hyou:20190406161454j:plain

こうなりました。たばこ・喫煙具専門小売業が約96%です。以下、自転車小売業、履物小売業(靴を除く)、骨とう品小売業、卵・鳥肉小売業、鮮魚小売業と続きました。

今回は以上です。

 次回は

 

www.crosshyou.info

 

です。

 

今回のR言語のコードです。


# head関数
head(df)

# summary関数
summary(df)

# 法人の割合
df4 <- df
df4$R_Corp <- df$Corp_Corp / df$Corp_Total
summary(df4$R_Corp)

# 個人の割合
df4$R_Indi <- df$Corp_Indi / df$Corp_Total
summary(df4$R_Indi)

# ヒストグラム
par(mfrow = c(1, 2))
hist(df4$R_Corp)
hist(df4$R_Indi)

# 年毎の法人割合の平均値
tapply(df4$R_Corp, df4$Year, mean)

# 年ごとの箱ひげ図
plot(factor(df4$Year), df4$R_Corp)

# 法人の割合 = a + b * 年
model1 <- lm(R_Corp ~ Year, data = df4)
summary(model1)

# 散布図と回帰直線
plot(df4$Year, df4$R_Corp)
abline(model1, col = "red", lty = 2)

# model1の残差プロット
plot(model1, which = 1)

# 一事業所当りの売上高
df4$Rev_Cor <- df4$Revenue / df4$Corp_Total
summary(df4$Rev_Cor)

# 法人の割合 = a + b1 * 年 + b2 * 一事業所当りの売上高
model2 <- lm(R_Corp ~ Year + Rev_Cor, data = df4)
summary(model2)

# 回帰残差のプロット
plot(model2, which = 1)

# model1とmodel2が有意に異なるか
anova(model1, model2)

# 一事業所当りの従業者数
df4$Stf_Cor <- df4$Staff / df4$Corp_Total
summary(df4$Stf_Cor)

# model3(model2 + 一事業所当りの従業者数)
model3 <- lm(R_Corp ~ Year + Rev_Cor + Stf_Cor, data = df4)
summary(model3)

# model4(model1 + 一事業所当りの従業者数)
model4 <- lm(R_Corp ~ Year + Stf_Cor, data = df4)
summary(model4)

# model4の回帰残差プロット
plot(model4, which = 1)

# 法人の割合の高い業種
head(rev(sort(tapply(df4$R_Corp, df4$Sector, mean))))

# 個人の割合の高い業種
head(rev(sort(tapply(df4$R_Indi, df4$Sector, mean))))