産業別就業人口比率と所得
/Users/fujimotokazuo/Dropbox/RStudio/e-statAPI/産業人口-所得.Rmd
kazuo.fujimoto2007@gmail.com
Fri Feb 24 11:10:22 JST 2023
1 概要
本レポートは、gacco「誰でも使えるオープンデータ」(誰でも使える統計オープンデータオフィシャルスタディノート: データサイエンス・オンライン講座 2021)の第１週３回目の講義のサンプルをAPIで取得しそれを回帰分析して動作を確認するスクリプトである。
データをAPIで取得する部分は、別スクリプト、e-statAPIフォルダの「Week1産業別-所得をAPIで.Rmd」によって取得された、「集計データ.rda」「集計データ2rda」を読み込んでいる。
集計データ2は、集計データにおいて、0/0の結果NaNとなった部分をゼロに置換している。これは、Excelの回帰分析処理が、そのように（非数値をゼロに置換）しているからで、RのlmでのNaNの処理とは異なっている。
library(tidyverse)
library(FactoMineR)
library(plotly)
library(htmlwidgets)
library(showtext)
showtext_auto(TRUE)2 データの読み込み
データは、e-statAPIをもちいて「Week1産業別-所得をAPIで.Rmd」で作成
元の処理で、NaNになっている部分をゼロに置き換えたもの。
3 C経済基盤とF労働をbind_rowsしてmerge
load("経済基盤.rda")
load("労働.rda")bind_rows(経済基盤 ,労働) %>% filter(調査年=="2015年度") %>%
pivot_wider(id_cols = 地域,names_from = 項目,values_from = value) -> bindCandF
bindCandF %>% names## [1] "地域" "C120110_課税対象所得"
## [3] "C120120_納税義務者数（所得割）" "F2201_第1次産業就業者数"
## [5] "F2211_第2次産業就業者数" "F2221_第3次産業就業者数"bindCandF %>% count(地域)## # A tibble: 1,741 × 2
## 地域 n
## <chr> <int>
## 1 三重県 いなべ市 1
## 2 三重県 亀山市 1
## 3 三重県 伊勢市 1
## 4 三重県 伊賀市 1
## 5 三重県 南伊勢町 1
## 6 三重県 名張市 1
## 7 三重県 四日市市 1
## 8 三重県 多気町 1
## 9 三重県 大台町 1
## 10 三重県 大紀町 1
## # … with 1,731 more rowsbindCandF %>% DT::datatable()bindCandF %>%
mutate(全産業従事者=F2201_第1次産業就業者数 + F2211_第2次産業就業者数 + F2221_第3次産業就業者数) %>%
mutate(１次産業割合=100*F2201_第1次産業就業者数/全産業従事者) %>%
mutate(２次産業割合=100*F2211_第2次産業就業者数/全産業従事者) %>%
mutate(３次産業割合=100*F2221_第3次産業就業者数/全産業従事者) %>%
mutate(一人あたり課税対象所得=1000*C120110_課税対象所得/`C120120_納税義務者数（所得割）`) -> 集計データ
集計データ %>% names## [1] "地域" "C120110_課税対象所得"
## [3] "C120120_納税義務者数（所得割）" "F2201_第1次産業就業者数"
## [5] "F2211_第2次産業就業者数" "F2221_第3次産業就業者数"
## [7] "全産業従事者" "１次産業割合"
## [9] "２次産業割合" "３次産業割合"
## [11] "一人あたり課税対象所得"集計データ %>% count(地域)## # A tibble: 1,741 × 2
## 地域 n
## <chr> <int>
## 1 三重県 いなべ市 1
## 2 三重県 亀山市 1
## 3 三重県 伊勢市 1
## 4 三重県 伊賀市 1
## 5 三重県 南伊勢町 1
## 6 三重県 名張市 1
## 7 三重県 四日市市 1
## 8 三重県 多気町 1
## 9 三重県 大台町 1
## 10 三重県 大紀町 1
## # … with 1,731 more rows#集計データ %>% head
#集計データ %>% tail(20)3.1 check
集計データ %>% select(地域) %>% unlist -> 地域API
length(地域API)# 1741## [1] 1741#length(地域Excel)# 1741集計データ %>% names## [1] "地域" "C120110_課税対象所得"
## [3] "C120120_納税義務者数（所得割）" "F2201_第1次産業就業者数"
## [5] "F2211_第2次産業就業者数" "F2221_第3次産業就業者数"
## [7] "全産業従事者" "１次産業割合"
## [9] "２次産業割合" "３次産業割合"
## [11] "一人あたり課税対象所得"write_excel_csv(集計データ,file = "CandFbyAPI.csv") # Excelでの計算確認用4 NaNをゼロに置き換え（Excelの回帰分析と合わせるため）
#集計データ %>% select(地域) %>% unlist() %>% str_detect(pattern="福島") -> Area_福島
#集計データ[Area_福島,]
集計データ %>% mutate(１次産業割合 = ifelse(is.nan(１次産業割合), 0, １次産業割合),
２次産業割合 = ifelse(is.nan(２次産業割合), 0, ２次産業割合),
３次産業割合 = ifelse(is.nan(３次産業割合), 0, ３次産業割合)) -> 集計データ2集計データ2 %>% count(地域)## # A tibble: 1,741 × 2
## 地域 n
## <chr> <int>
## 1 三重県 いなべ市 1
## 2 三重県 亀山市 1
## 3 三重県 伊勢市 1
## 4 三重県 伊賀市 1
## 5 三重県 南伊勢町 1
## 6 三重県 名張市 1
## 7 三重県 四日市市 1
## 8 三重県 多気町 1
## 9 三重県 大台町 1
## 10 三重県 大紀町 1
## # … with 1,731 more rows5 データ（.rda）のsave
#save(集計データ,file="集計データ.rda") # NaN をそのまま
#save(集計データ2,file="集計データ2.rda")# NaNをゼロに置換#load("集計データ.rda")
#load("集計データ2.rda")5.1 データの概要を確認する
5.1.1 NaN をゼロ置換してない元のデータ
集計データ %>% summary()## 地域 C120110_課税対象所得 C120120_納税義務者数（所得割）
## Length:1741 Min. :2.824e+05 Min. : 100
## Class :character 1st Qu.:8.155e+06 1st Qu.: 3212
## Mode :character Median :2.675e+07 Median : 10082
## Mean :1.055e+08 Mean : 32095
## 3rd Qu.:7.797e+07 3rd Qu.: 27072
## Max. :7.084e+09 Max. :1777500
##
## F2201_第1次産業就業者数 F2211_第2次産業就業者数 F2221_第3次産業就業者数
## Min. : 0 Min. : 0 Min. : 0
## 1st Qu.: 357 1st Qu.: 984 1st Qu.: 2313
## Median : 763 Median : 2971 Median : 7256
## Mean : 1276 Mean : 7996 Mean : 22754
## 3rd Qu.: 1622 3rd Qu.: 8276 3rd Qu.: 19150
## Max. :15563 Max. :324156 Max. :1233147
##
## 全産業従事者 １次産業割合 ２次産業割合 ３次産業割合
## Min. : 0 Min. : 0.000 Min. : 1.539 Min. :19.80
## 1st Qu.: 3934 1st Qu.: 2.975 1st Qu.:19.752 1st Qu.:56.03
## Median : 11478 Median : 8.051 Median :25.640 Median :62.21
## Mean : 32026 Mean :11.121 Mean :25.814 Mean :63.06
## 3rd Qu.: 29440 3rd Qu.:16.361 3rd Qu.:31.529 3rd Qu.:70.36
## Max. :1565064 Max. :77.061 Max. :69.892 Max. :93.40
## NA's :5 NA's :5 NA's :5
## 一人あたり課税対象所得
## Min. : 1984254
## 1st Qu.: 2491690
## Median : 2710458
## Mean : 2801271
## 3rd Qu.: 2990759
## Max. :10232188
## 5.1.2 NaN をゼロ置換したデータ
集計データ2 %>% summary()## 地域 C120110_課税対象所得 C120120_納税義務者数（所得割）
## Length:1741 Min. :2.824e+05 Min. : 100
## Class :character 1st Qu.:8.155e+06 1st Qu.: 3212
## Mode :character Median :2.675e+07 Median : 10082
## Mean :1.055e+08 Mean : 32095
## 3rd Qu.:7.797e+07 3rd Qu.: 27072
## Max. :7.084e+09 Max. :1777500
## F2201_第1次産業就業者数 F2211_第2次産業就業者数 F2221_第3次産業就業者数
## Min. : 0 Min. : 0 Min. : 0
## 1st Qu.: 357 1st Qu.: 984 1st Qu.: 2313
## Median : 763 Median : 2971 Median : 7256
## Mean : 1276 Mean : 7996 Mean : 22754
## 3rd Qu.: 1622 3rd Qu.: 8276 3rd Qu.: 19150
## Max. :15563 Max. :324156 Max. :1233147
## 全産業従事者 １次産業割合 ２次産業割合 ３次産業割合
## Min. : 0 Min. : 0.000 Min. : 0.00 Min. : 0.00
## 1st Qu.: 3934 1st Qu.: 2.950 1st Qu.:19.71 1st Qu.:55.97
## Median : 11478 Median : 8.033 Median :25.61 Median :62.19
## Mean : 32026 Mean :11.089 Mean :25.74 Mean :62.88
## 3rd Qu.: 29440 3rd Qu.:16.296 3rd Qu.:31.51 3rd Qu.:70.35
## Max. :1565064 Max. :77.061 Max. :69.89 Max. :93.40
## 一人あたり課税対象所得
## Min. : 1984254
## 1st Qu.: 2491690
## Median : 2710458
## Mean : 2801271
## 3rd Qu.: 2990759
## Max. :102321885.1.3 値のチェック！
集計データ2[集計データ2$一人あたり課税対象所得==0,]## # A tibble: 0 × 11
## # … with 11 variables: 地域 <chr>, C120110_課税対象所得 <dbl>,
## # C120120_納税義務者数（所得割） <dbl>, F2201_第1次産業就業者数 <dbl>,
## # F2211_第2次産業就業者数 <dbl>, F2221_第3次産業就業者数 <dbl>,
## # 全産業従事者 <dbl>, １次産業割合 <dbl>, ２次産業割合 <dbl>,
## # ３次産業割合 <dbl>, 一人あたり課税対象所得 <dbl>集計データ2 %>% select(地域,１次産業割合,一人あたり課税対象所得) %>% DT::datatable()5.2 ヒストグラムで分布を確認する
集計データ2 %>% select(1,8:10) %>% pivot_longer(cols = -`地域`,names_to = "産業割合",values_to = "割合") -> 産業割合.long
産業割合.long %>% ggplot(aes(x=割合,fill=産業割合)) + geom_histogram(position="identity",alpha=0.5,binwidth = 1) -> p
p <- p + ggtitle("人口に対する各産業区分分布（2015社会・人口体系）") # position="identity"を指定しないと重ねた度数のヒストグラム
p
5.3 EDAtoolのexploreでの集計は参考まで
6 変数の分布を確認する
explore::report(集計データ2,output_dir = ".",output_file = "report.html")explore::explore(集計データ2)7 回帰分析
7.1 グラフ化
7.2 第１次産業比率と所得
集計データ2 %>% ggplot(aes(x=１次産業割合,y=一人あたり課税対象所得,label=地域)) +geom_point() + geom_smooth(method=lm) +
ggtitle("第１次産業比率と所得") -> p
ggplotly(p) -> pp## `geom_smooth()` using formula 'y ~ x'pp#saveWidget(pp,"graph1.html")7.2.1 回帰式を求める
res.lm <- lm(一人あたり課税対象所得~１次産業割合,data = 集計データ2)
summary(res.lm) -> res.sum.lm
res.sum.lm##
## Call:
## lm(formula = 一人あたり課税対象所得 ~ １次産業割合,
## data = 集計データ2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -971023 -257670 -96864 136169 7211743
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3022236 16990 177.89 <2e-16 ***
## １次産業割合 -19926 1121 -17.78 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 483200 on 1739 degrees of freedom
## Multiple R-squared: 0.1538, Adjusted R-squared: 0.1533
## F-statistic: 316 on 1 and 1739 DF, p-value: < 2.2e-16res.sum.lm$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3022236.4 16989.593 177.88751 0.000000e+00
## １次産業割合 -19925.8 1120.931 -17.77612 4.358081e-657.2.2 NaN未処理のデータにlm
res.lm <- lm(一人あたり課税対象所得~１次産業割合, 集計データ)
summary(res.lm)##
## Call:
## lm(formula = 一人あたり課税対象所得 ~ １次産業割合,
## data = 集計データ)
##
## Residuals:
## Min 1Q Median 3Q Max
## -968956 -257772 -97702 135716 7214033
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3019936 17041 177.21 <2e-16 ***
## １次産業割合 -19815 1123 -17.65 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 483200 on 1734 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.1523, Adjusted R-squared: 0.1518
## F-statistic: 311.5 on 1 and 1734 DF, p-value: < 2.2e-16summary(res.lm) -> res.sum.lm
res.sum.lm$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3019936.04 17041.247 177.21332 0.000000e+00
## １次産業割合 -19814.76 1122.723 -17.64884 3.058378e-647.3 第二次産業比率と所得
集計データ2 %>% ggplot(aes(x=２次産業割合,y=一人あたり課税対象所得,label=地域)) +geom_point() + geom_smooth(method=lm) +
ggtitle("第２次産業比率と所得") -> p
ggplotly(p)## `geom_smooth()` using formula 'y ~ x'7.3.1 回帰式を求める
res.lm <- lm(一人あたり課税対象所得~２次産業割合, 集計データ2)
summary(res.lm) -> res.sum.lm
res.sum.lm##
## Call:
## lm(formula = 一人あたり課税対象所得 ~ ２次産業割合,
## data = 集計データ2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -825664 -306050 -90902 188054 7398818
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2858334 40652 70.312 <2e-16 ***
## ２次産業割合 -2217 1502 -1.476 0.14
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 525000 on 1739 degrees of freedom
## Multiple R-squared: 0.001251, Adjusted R-squared: 0.0006771
## F-statistic: 2.179 on 1 and 1739 DF, p-value: 0.1401res.sum.lm$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2858333.923 40652.423 70.311526 0.0000000
## ２次産業割合 -2216.883 1501.799 -1.476152 0.14008427.4 第３次産業比率と所得
集計データ2 %>% ggplot(aes(x=３次産業割合,y=一人あたり課税対象所得,label=地域)) +geom_point() + geom_smooth(method=lm) +
ggtitle("第３次産業比率と所得") -> p
ggplotly(p)## `geom_smooth()` using formula 'y ~ x'7.4.1 回帰式を求める
res.lm <- lm(一人あたり課税対象所得~３次産業割合, 集計データ2)
summary(res.lm) -> res.sum.lm
res.sum.lm##
## Call:
## lm(formula = 一人あたり課税対象所得 ~ ３次産業割合,
## data = 集計データ2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1037526 -274242 -76065 178974 6948926
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1624939 69921 23.24 <2e-16 ***
## ３次産業割合 18707 1096 17.06 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 486200 on 1739 degrees of freedom
## Multiple R-squared: 0.1434, Adjusted R-squared: 0.1429
## F-statistic: 291.1 on 1 and 1739 DF, p-value: < 2.2e-16res.sum.lm$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1624939.31 69920.839 23.23970 2.894111e-104
## ３次産業割合 18706.57 1096.365 17.06236 1.784700e-608 重回帰分析
8.1 NaNをゼロ処理したデータ
res.lm <- lm(一人あたり課税対象所得~１次産業割合+２次産業割合+３次産業割合, 集計データ2)
summary(res.lm)##
## Call:
## lm(formula = 一人あたり課税対象所得 ~ １次産業割合 +
## ２次産業割合 + ３次産業割合, data = 集計データ2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1029851 -244616 -68160 154409 6938774
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3392143 210151 16.141 < 2e-16 ***
## １次産業割合 -24246 2319 -10.457 < 2e-16 ***
## ２次産業割合 -14281 2360 -6.051 1.76e-09 ***
## ３次産業割合 725 2155 0.336 0.737
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 469900 on 1737 degrees of freedom
## Multiple R-squared: 0.2007, Adjusted R-squared: 0.1993
## F-statistic: 145.4 on 3 and 1737 DF, p-value: < 2.2e-168.2 NaNをゼロ処理しないデータ
- ３次産業割合のNAはなぜ発生したのか…
res.lm <- lm(一人あたり課税対象所得~１次産業割合+２次産業割合+３次産業割合, 集計データ)
summary(res.lm)##
## Call:
## lm(formula = 一人あたり課税対象所得 ~ １次産業割合 +
## ２次産業割合 + ３次産業割合, data = 集計データ)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1029851 -245149 -68168 154515 6938774
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3464638 47706 72.625 <2e-16 ***
## １次産業割合 -24970 1209 -20.651 <2e-16 ***
## ２次産業割合 -15006 1509 -9.942 <2e-16 ***
## ３次産業割合 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 470100 on 1733 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.198, Adjusted R-squared: 0.1971
## F-statistic: 213.9 on 2 and 1733 DF, p-value: < 2.2e-169 主成分分析（PCA）
res.PCA <- PCA(集計データ[,8:10])## Warning in PCA(集計データ[, 8:10]): Missing values are imputed by the mean of
## the variable: you should use the imputePCA function of the missMDA package
res.PCA2 <- PCA(集計データ2[,8:10])
参考文献
誰でも使える統計オープンデータオフィシャルスタディノート:
データサイエンス・オンライン講座. 2021. 改訂第2版 ed. Tokyo:
日本統計協会.