#library(xlsx)# R3.4.0 でrjavaが動かない？
library(readxl)# Excel の読み込みにこちらを使う(^^;)
library(FactoMineR)
library(factoextra)
library(ggpubr)
library(tidyverse)
library(stringr)# https://heavywatal.github.io/rstats/stringr.html
#library(ggbiplot)このデータ（刑法犯認知件数）をつかって主成分分析を行うのは、Tokyo.R#60（https://atnd.org/events/87238） での、@kilometer00 さんのプレゼンのDEMOで教えていただきました。
これを事例にすると、合成得点をだす「主成分分析」と組み合わせパターンの関係（構造）をみる「対応分析」の特徴がわかりやすいかと。
データを警察白書H28から取得
https://www.npa.go.jp/hakusyo/h28/data.html
#.d0 <- read.xlsx(file = "01.xls",1,encoding="cp932",startRow =7,stringsAsFactors=FALSE)
.d0 <- read_excel("01.xls",skip=6)
colnames(.d0)[1:2] <- c("地方区分","都道府県")
names(.d0)## [1] "地方区分" "都道府県" "認知件数" "検挙件数" "検挙人員"
## [6] "認知件数__1" "検挙件数__1" "検挙人員__1" "認知件数__2" "検挙件数__2"
## [11] "検挙人員__2" "認知件数__3" "検挙件数__3" "検挙人員__3" "認知件数__4"
## [16] "認知件数__5" "検挙件数__4" "検挙人員__4" "認知件数__6" "検挙件数__5"
## [21] "検挙人員__5" "認知件数__7" "検挙件数__6" "検挙人員__6".d0[15,2] <- "東京" # 東京だけセルがずれている
#.d0 %>% filter(都道府県!="計")
#.d0t <- read.xlsx(file = "01.xls",1,encoding="cp932",startRow = 4,endRow = 5,stringsAsFactors=FALSE)
# 犯罪項目見出しを取得
.d0t <- read_excel("01.xls",skip=3,n_max=5,col_types = rep("text",24))
pos <- c(seq(3,9,3),c(13,16,19,22))
c_label <- .d0t[1,pos]
.d0 %>% filter(都道府県!="計") %>% select(c(2,pos)) -> 認知件数
colnames(認知件数) <- c("都道府県",names(unlist(c_label)))
colnames(認知件数)## [1] "都道府県" "刑法犯総数" "凶悪犯" "粗暴犯"
## [5] "窃盗犯" "知能犯" "風俗犯" "その他の刑法犯"rownames(認知件数) <- unlist(認知件数[,1])#PCAのためにrownamesをつける## Warning: Setting row names on a tibble is deprecated.colnames(認知件数)[8] <- "その他"# 「その他の刑法犯」を「その他」に略記主成分分析１ 認知件数
#res.pca <- princomp(認知件数[,2:ncol(認知件数)])
#biplot(res.pca)
#ggbiplot(res.pca) + theme_gray(base_family = "sans")
.d <- data.frame(認知件数)
res.PCA <- PCA(.d[,2:ncol(認知件数)],graph=FALSE)
fviz_pca_biplot(res.PCA,font.family = "sans",repel =FALSE,title="都道府県別刑法犯認知件数/人口比 主成分分析") + theme(text = element_text(family = "sans"))
主成分分析２ 人口比率で検討する
認知件数[,1:ncol(認知件数)]## # A tibble: 51 × 8
## 都道府県 刑法犯総数 凶悪犯 粗暴犯 窃盗犯 知能犯 風俗犯 その他
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 札幌 24749 112 1379 4697 742 750 4863
## 2 函館 2407 23 134 606 91 49 528
## 3 旭川 3497 11 299 1054 112 59 753
## 4 釧路 3641 28 266 837 126 91 822
## 5 北見 1163 10 95 304 55 9 266
## 6 青森 5486 35 421 1536 303 90 1028
## 7 岩手 4884 23 264 1518 244 81 851
## 8 宮城 17742 78 922 3785 989 229 2810
## 9 秋田 3154 24 303 1231 220 35 382
## 10 山形 5014 19 649 1886 233 44 641
## # ... with 41 more rows#これに人口データをつけるただ、これは認知件数をそのまま評価しているので、人口の大小の影響をみているだけでもあるので、人口比でみることにする。次へ。
人口データを読み込む
警察庁のデータとあわせるために、H27年度の２を使う http://www.e-stat.go.jp/SG1/estat/List.do?lid=000001163203
北海道が、なぜか市部の統計なので、その部分の人口もとってくる。それにしても、#小樽、室蘭はない？？ http://www.e-stat.go.jp/SG1/estat/GL08020103.do?_toGL08020103_&tclassID=000001077438&cycleCode=0&requestSender=estat
#.dj0 <- read.xlsx(file = "002.xls",1,encoding="cp932",startRow =18,stringsAsFactors=FALSE)
.dj0 <- read_excel("002.xls",skip=17,col_types = c(rep("text",24)))
.dj1 <- .dj0[1:47,c(3,6)]
colnames(.dj1) <- c("都道府県","人口")
.dj2 <- data.frame(.dj1)
.dj2[,2] <- as.numeric(.dj2[,2])*1000
.dj2[,1] <- str_replace(.dj2[,1],"県","")
.dj2[,1] <- str_replace(.dj2[,1],"東京都","東京")
.dj2[,1] <- str_replace(.dj2[,1],"大阪府","大阪")
.dj2[,1] <- str_replace(.dj2[,1],"京都府","京都")
#.dj2[,1] <- str_replace(.dj2[,1],"京","京都")
hkd_pop <- data.frame(都道府県 = c("札幌","函館","旭川","釧路","北見"),
人口 = #c(1952,266,340,175,121))
c(1952356,265979,339605,174742,121226)
)#小樽、室蘭はない？？
rownames(hkd_pop) <- hkd_pop[,1]
人口 <- rbind(hkd_pop,.dj2)
人口[,1] <- stringr::str_trim(人口[,1],side="both")# 都道府県カラムにある空白を削除
認知件数 %>% left_join(人口) -> 認知件数2#人口 %>% right_join(認知件数)## Joining, by = "都道府県"# > mutate(mtcars, v1 = hp/cyl, v2 = v1/wt)
認知件数2 %>% dplyr::transmute(都道府県,
凶悪=凶悪犯/人口,
粗暴=粗暴犯/人口,
窃盗=窃盗犯/人口,
知能=知能犯/人口,
風俗=風俗犯/人口,
他=その他/人口) -> 認知件数p
.d <- data.frame(認知件数p)
rownames(.d) <- unlist(認知件数p[,1])
res.PCA <- PCA(.d[,2:ncol(認知件数p)],graph=FALSE)## Warning in PCA(.d[, 2:ncol(認知件数p)], graph = FALSE): Missing values are
## imputed by the mean of the variable: you should use the imputePCA function
## of the missMDA packagefviz_pca_biplot(res.PCA,font.family = "sans",repel =FALSE)
対応分析
認知件数 %>% select(3:8) -> .d2
rownames(.d2) <- unlist(認知件数[,1])## Warning: Setting row names on a tibble is deprecated.res.CA <- CA(.d2,graph = FALSE)
fv <- fviz_ca_biplot(res.CA,arrows = c(FALSE,TRUE),map="colgreen",font.family = "sans",title="都道府県別刑法犯認知件数 対応分析")
fv + theme(text = element_text(family = "sans"))
因子分析
res.FA <- factanal(.d2,factors=3)
res.FA##
## Call:
## factanal(x = .d2, factors = 3)
##
## Uniquenesses:
## 凶悪犯 粗暴犯 窃盗犯 知能犯 風俗犯 その他
## 0.005 0.022 0.005 0.009 0.082 0.045
##
## Loadings:
## Factor1 Factor2 Factor3
## 凶悪犯 0.745 0.437 0.499
## 粗暴犯 0.553 0.620 0.537
## 窃盗犯 0.471 0.784 0.397
## 知能犯 0.527 0.545 0.645
## 風俗犯 0.752 0.478 0.352
## その他 0.572 0.561 0.560
##
## Factor1 Factor2 Factor3
## SS loadings 2.253 2.030 1.549
## Proportion Var 0.376 0.338 0.258
## Cumulative Var 0.376 0.714 0.972
##
## The degrees of freedom for the model is 0 and the fit was 0.0042barplot(res.FA$loadings[,1:2],beside = TRUE)
plot(res.FA$loadings[,1:2],type="n")
text(res.FA$loadings[,1:2],colnames(.d2)) 