manuscript01
Yutaka Moteki
8/27/2020
1. Rによる二次的データの分析
1.1 はじめに
Rの使い方をまとめてみようと思う。Rの使い方についてはおよそ6カ月間毎日取り組んだので、自分が使う可能性がある手法についてはだいたいマスターすることができたはずである。理解したことをまとめてみて検討してみようと思う。Rと統計学に関する知識を整理してみる。
また、このファイル（Rmd形式）を編集することで、RStudioが原稿などを書くのに使えるかどうかを確認してみたい。グラフや表のある原稿はRStudioが便利のようである。
（html_notebookを使って見たのだが使い方がややこしいので、元のhtml_documentに戻した。）
（次の行は空けておかないと##が有効にならない。空白行は改行と同じ効果。改行は半角スペース２個でも可能。）
1.2 long formatの活用
long formatは、縦持ち形式とか縦長形式とか訳される。これについて最初に検討してみよう。 いままで使ってきて、視覚的にとらえやすいのは変数列が横長に広がるwide format1であるが、long formatの方が便利な面があるということなので試している。wide formatをlong formatに変換する関数について確認してみる（gatherとpivot_longer）。（注は、ドキュメントの最後に表示される。）
ここで、 RStudioで文章を書いたり数式を示したりする場合の約束事について確認してみる。 まず、数式で使う記号の確認。**が^として使える。Rでの実行結果で確認できるように、２の3乗は次のどちらでもよい。（この部分は、実際には次のように入力している。）
\*\*が^として使える。**２の3乗**は次のどちらでもよい。
> 2**3
[1] 8
> 2^3
[1] 8字体について。 2つのアスタリスク(*)で囲むとボールドになる。ひとつのアスタリスクで囲むとイタリック。
記号**を使うとボールド体の指定のはじまりと解釈されてしまうばあいがある（後から**が出てくる場合）。それを避けるためにはその前に記号\を前に1つ入れるとよい。（この文の実際の入力は以下のとおり。）
記号\*\*を使うとボールド体の指定のはじまりと解釈されてしまうばあいがある（後から\*\*が出てくる場合）。それを避けるためにはその前に記号\\を前に*1つ*入れるとよい。「`」の使い方。この記号で囲むと囲まれた部分をコマンドなどとして表記できる。複数行にわたるものの場合は、最初と最後に、この記号を３個並べた「```」を使う。
1.3 gather関数の使い方
wide formatをlong formatに変換するこの関数は、次のような使い方をする。
gather(df,"key","value",-x)
gather(df,"key","value",y1,y2,y3)
gather(
data,
key = "key",
value = "value",
...,
na.rm = FALSE,
convert = FALSE,
factor_key = FALSE
)gather関数を使うとwide formatをlong formatに変換することができる。これはtydyrパッケージに含まれている。 (1)の”key”と”value”のところは、別の名前をつけてよい。新変数であり、日本語を使ってもよい。新変数は引用符で囲むことが推奨されている。
(1)の場合、yがこの形式では、設定されていない。マイナス記号を付けた変数は、インデックス変数として扱われ、それ以外が新変数の”key”と”value”を使って集められる。
(2)の場合、y1,y2,y3以外は、インデックス変数として扱われる。 y1,y2,y3をy1:y3と表現できる。どちらで表記してもよい。ここに出てこない変数(ここではx)はインデックス変数として扱われる。
想定されているのは以下のようなデータフレーム（オブジェクト名:df）である。これはx,y1,y2,y3の3変数からなるデータフレームである。
library(tidyr)
library(ggplot2)
x <- 1:10
y1 <- x**1.5
y2 <- x**2
y3 <- x**2.5
df <- data.frame(x,y1,y2,y3)
knitr::kable(df)| x | y1 | y2 | y3 |
|---|---|---|---|
| 1 | 1.000000 | 1 | 1.000000 |
| 2 | 2.828427 | 4 | 5.656854 |
| 3 | 5.196152 | 9 | 15.588457 |
| 4 | 8.000000 | 16 | 32.000000 |
| 5 | 11.180340 | 25 | 55.901699 |
| 6 | 14.696938 | 36 | 88.181631 |
| 7 | 18.520259 | 49 | 129.641814 |
| 8 | 22.627417 | 64 | 181.019336 |
| 9 | 27.000000 | 81 | 243.000000 |
| 10 | 31.622777 | 100 | 316.227766 |
gather関数を使い横長のdfを縦長のdf_longに変えるには以下のようにする。変換されたデータフレームはdf_longという名前にしている。長くなるので、13行目まで表示する。
df_long <- gather(df,"key","value",-x)
knitr::kable(df_long[1:13,])| x | key | value |
|---|---|---|
| 1 | y1 | 1.000000 |
| 2 | y1 | 2.828427 |
| 3 | y1 | 5.196152 |
| 4 | y1 | 8.000000 |
| 5 | y1 | 11.180340 |
| 6 | y1 | 14.696938 |
| 7 | y1 | 18.520259 |
| 8 | y1 | 22.627417 |
| 9 | y1 | 27.000000 |
| 10 | y1 | 31.622777 |
| 1 | y2 | 1.000000 |
| 2 | y2 | 4.000000 |
| 3 | y2 | 9.000000 |
あるいは、2番目の形式で以下のようにしても同じ結果である。
df_long <- gather(df,"key","value",y1,y2,y3)
#knitr::kable(df_long[1:13,])調べてみてわかったことだが、
gather関数はすでに開発終了(retired
lifecyle)で、以下のように、pivot-longerに移行することが推奨されている。
Gather columns into key-value pairs
Description
Retired lifecycle
Development on gather() is complete, and for new code we recommend switching to pivot_longer(), which is easier to use, more featureful, and still under active development. df %>% gather("key", "value", x, y, z) is equivalent to df %>% pivot_longer(c(x, y, z), names_to = "key", values_to = "value")
See more details in vignette("pivot").
Usage
gather(
data,
key = "key",
value = "value",
...,
na.rm = FALSE,
convert = FALSE,
factor_key = FALSE
)次のような式が推奨されている。これに従って書き換えてみると同じようなが得られる（並べ方には違いがある）。gather関数の場合と同様に、“key”と”value”のところは変更可能であり、日本語で”凡例”,“値”などとしてもよい。(keyのところは凡例にすれば凡例が凡例と表示されるので便利であるが、複数の変数をまとめる新変数に「凡例」という名前が付くのは感覚的には受け入れにくいことである。その意味ではあとで、labsのなかで指定する方がよいかもしれない。)
df %>% pivot_longer(c(x, y, z), names_to = "key", values_to = "value")library(tidyr)
df_long <-df %>%
pivot_longer(c(y1, y2, y3), names_to = "key", values_to = "value")
df_long## # A tibble: 30 × 3
## x key value
## <int> <chr> <dbl>
## 1 1 y1 1
## 2 1 y2 1
## 3 1 y3 1
## 4 2 y1 2.83
## 5 2 y2 4
## 6 2 y3 5.66
## 7 3 y1 5.20
## 8 3 y2 9
## 9 3 y3 15.6
## 10 4 y1 8
## # ℹ 20 more rowspivot_longerを使った場合には、xを軸にしてy1,y2,y3の値が示されている。
データを用いてグラフを描く場合、 wide
formatの新変数名に対応してggplot2でのaesの設定を変える。
図１では、ggplot2のaesでx=x,y=value,col=keyと設定している。
x軸、y軸、凡例のタイトルについてlabsの中で設定してはいない。
図２では、凡例及び値という変数をpivot_longer関数の中で設定し、それらを使ってggplot2のaesの中でx=x,y=値,col=凡例と指定している。y軸のラベル、凡例のタイトルにそれらが自動的に反映している。
library(patchwork)
p1 <- ggplot(df_long,aes(x=x,y=value,col=key))+
geom_line()+
labs(title="図1",subtitle = "aes(x=x,y=value,col=key)")+
theme_minimal(base_family="HiraKakuProN-W3")
df_long <-df %>%
pivot_longer(c(y1, y2, y3), names_to = "凡例", values_to = "値")
p2 <- ggplot(df_long,aes(x=x,y=値,col=凡例))+
geom_line()+
labs(title="図2",subtitle = "aes(x=x,y=値,col=凡例)")+
theme_minimal(base_family="HiraKakuProN-W3")df_long <- gather(df,"凡例","y",y1:y3)
p3 <- ggplot(df_long,aes(x,y,linetype=凡例))+
geom_line()+
labs(title="図3",subtitle = "aes(x,y,linetype=凡例)")+
theme_minimal(base_family="HiraKakuProN-W3")
p4 <- ggplot(df_long,aes(x,y,linetype=凡例))+
geom_line(col="blue")+
labs(title="図4",subtitle = "aes(x,y,linetype=凡例)")+
theme_minimal(base_family="HiraKakuProN-W3")
図３では、新変数を凡例、yと名付け、gather(df,"凡例","y",y1:y3)を実行し、さらにggplot2のaesでaes(x,y,linetype=凡例)と指定している。これは、aes(x=x,y=y,linetype=凡例)ということと同じである。この図3は、モノクロで線の種類を使い分けている。図４では、さらに、aesという形式ではなく、geom_lineについて、geom_line(col="blue")として青色を指定している。１色だけであるが、様々な色を線に指定することができる。
1.4 ggplot2でグラフを描く
ggplot2を使うと、さまざまなグラフを自由に描くことできる。方程式を使って曲線と直線を同じグラフに描いてみたり、軸を反転させたりしてみよう。(図2は図1の軸を反転させたものである。)
図１では、新変数としてbucketとyを使っている。
そしてggplot2のaesでx,y,col=bucketと指定している。
しかし、labsの中でcol="凡例"と設定しているので凡例は凡例と表示される。
labsでタイトル、軸、凡例などの名前を指定できる。
y軸のラベルはここでは設定されていないが、その上のaesの設定でy=yということなのでy軸にはｙが表示される。
図2は、coord_flip()を使ってｘ軸とｙ軸を交換している。また、モノクロで線を描いている。新変数は凡例とyであり、linetype=凡例とaesの中で設定している。なお、図2では、captionを使い、グラフの右下に説明を付け加えている。キャプションは以下のように設定している。2行に分けて表示されるように「\n」で改行をおこなっている。
caption = "図1の軸を交換したもの。linetypeに凡例\nをわりあて、点ではなく線を描く。"x <- ppoints(100)
y1 <- x**(-0.5)
y2 <- -8*x**2+14.23
y3 <- -14*x + 14.23
df <- data.frame(x,y1,y2,y3)
df_long <- gather(df,"bucket","y",-x)
library(patchwork)
p1 <- ggplot(df_long,aes(x,y,col=bucket))+
geom_point()+
labs(title="図1",subtitle="y1 = x**(-0.5)\ny2 = -8*x**2+14.23\ny3 = -14*x + 14.23",col="凡例")+
theme_minimal(base_family="HiraKakuProN-W3")
df_long <- gather(df,"凡例","y",-x) # 新変数を"凡例","ｙ"と設定。
p2 <- ggplot(df_long,aes(x,y,linetype=凡例))+
geom_line()+
labs(title="図2",subtitle="y1 = x**(-0.5)\ny2 = -8*x**2+14.23\ny3 = -14*x + 14.23",caption = "図1の軸を交換したもの。linetypeに凡例\nをわりあて、点ではなく線を描く。")+
theme_minimal(base_family="HiraKakuProN-W3")+
coord_flip()
1.5 SSDSE
colnamesで元のデータフレームの変数名を変更すればグラフの凡例にも反映する(図1、図2)。
df <- data.frame(pop,d22,OutOverIn)
colnames(df) <- c("pop","財政力指数","転出超過率")
df_long <- gather(df,"bucket","y",-pop)以下のcolour(colでもよい)の指定に注意。
labs(title = "人口と財政力指数、転出超過率",subtitle = "2つのデータ系列の同時表示",x="人口",y="値",colour="凡例")+labsの中のこの場合は、凡例のタイトルが指定される。ここでは「凡例」としたのでそのように表示されている。指定しなければ新変数bucketとなる。
図３では、labsの中でcaptionを指定している。その内容はグラフの右下に表示される。
library(tidyr)
library(ggplot2)
dat <- read.csv("./data/dat_ssdse.csv")
pop <- dat$A1101
d22 <- dat$D2201
OutOverIn <- dat$OutOverIn
df <- data.frame(pop,d22,OutOverIn)
colnames(df) <- c("pop","財政力指数","転出超過率")
df_long <- gather(df,"bucket","y",-pop)
library(patchwork)
p1 <- ggplot(df_long,aes(pop,y,col=bucket))+
geom_point()+
labs(title = "図１. 人口と財政力指数、転出超過率",subtitle = "2つのデータ系列の同時表示",x="人口",y="値",colour="凡例")+
theme_minimal(base_family="HiraKakuProN-W3")
p2 <- ggplot(df_long,aes(pop,y,col=bucket))+
geom_point()+
xlim(0,1e+06)+
labs(title = "図２. 人口と財政力指数、転出超過率",subtitle = "2つのデータ系列の同時表示",x="人口",y="値",col="凡例")+
theme_minimal(base_family="HiraKakuProN-W3")
df <- data.frame(pop,d22)
df_long <- gather(df,"key","value",-pop)
p3 <- ggplot(df_long,aes(x=pop,y=value,colour="財政力指数"))+
geom_point(colour="blue")+
xlim(0,1e+06)+
labs(title = "図３. 人口と財政力指数",subtitle = "図２の赤い部分だけを表示",x="人口",y="財政力指数",colour="凡例",caption = "2015年度のデータ")+
theme_minimal(base_family="HiraKakuProN-W3")
library(tidyr)
library(ggplot2)
dat <- read.csv("./data/dat_ssdse.csv")
x <- dat$Dnsty # 可住地人口密度
y1 <- dat$D2201 # 財政力指数
y2 <- dat$r_OnePerson # 単独世帯率
y3 <- dat$OutOverIn # 転出超過率
df <- data.frame(x,y1,y2,y3)
df_long <- gather(df,"凡例","y",-x)
ggplot(df_long,aes(x,y,col=凡例))+
geom_point()+
labs(title = "可住地人口密度と財政力指数等",subtitle = "y1: 財政力指数\ny2: 単独世帯率\ny3: 転出超過率",x="可住地人口密度",y="y")+
theme_minimal(base_family="HiraKakuProN-W3")
2. 散布図を描く
2.1 散布図(全国)
点の色は、geom_point(colour="red")などとして指定できる。
library(tidyr)
library(ggplot2)
library(ggrepel)
dat <- read.csv("./data/dat_ssdse.csv")
x <- dat$Dnsty # 可住地人口密度
y1 <- dat$D2201 # 財政力指数
df <- data.frame(x,y1)
df_long <- gather(df,"凡例","y",-x)
ggplot(df_long,aes(x,y))+
geom_point(colour="seagreen")+
labs(title = "全国",subtitle = "y: 財政力指数",x="x: 可住地人口密度",y="y")+
theme_minimal(base_family="HiraKakuProN-W3")
library(tidyr)
library(ggplot2)
dat <- read.csv("./data/dat_ssdse.csv")
x <- dat$Dnsty # 可住地人口密度
y1 <- dat$D2201 # 財政力指数
df <- data.frame(x,y1)
df_long <- gather(df,"凡例","y",-x)
ggplot(df_long,aes(x,y,linetype="凡例"))+
geom_point()+
xlim(0,10)+
labs(title = "全国（xの範囲を限定）",subtitle = "y: 財政力指数",x="x: 可住地人口密度",y="y")+
theme_minimal(base_family="HiraKakuProN-W3")
2.2 雛形（全国用）
以下は、colnames(dat)を実行した結果であり、変数名のリストである。これによって各変数がデータフレームの何列目にあるかを調べる。
## [1] "X" "Code" "Prefecture"
## [4] "Municipality" "A1101" "A6108"
## [7] "B1103" "A5101" "A5102"
## [10] "D2201" "F1102" "F2701"
## [13] "A710201" "A710202" "A1303"
## [16] "A4101" "A4200" "A710101"
## [19] "A810101" "A810102" "A810103"
## [22] "A810104" "A810105" "D2205"
## [25] "D2206" "Dnsty" "In"
## [28] "Out" "OutOverIn" "sc"
## [31] "r_instHm" "r_nuclear" "r_nuclear_relativeH"
## [34] "r_OnePerson" "aged" "naturalInc"csvファイルを読み込み、datというオブジェクトとして保存する。Codeとともに、2変数を取り出す。gather式の中で、Codeとx(横軸となる変数)をインデックス変数として指定する。以下の例では、新変数「凡例」に 「sc」が入れられる。
library(tidyr)
library(ggplot2)
dat <- read.csv("./data/dat_ssdse.csv")
colnames(dat)[c(2,26,30)]## [1] "Code" "Dnsty" "sc"df_long <- gather(dat[,c(2,26,30)],"凡例","y",-Code,-Dnsty)
head(df_long)## Code Dnsty 凡例 y
## 1 1100 44.474828 sc 44.14595
## 2 1202 21.566448 sc 87.57994
## 3 1203 15.223374 sc 80.90107
## 4 1204 9.649514 sc 91.29442
## 5 1205 20.467761 sc 85.60854
## 6 1206 5.653251 sc 87.12453tail(df_long)## Code Dnsty 凡例 y
## 1736 47360 1.3077586 sc 99.32341
## 1737 47361 1.9431220 sc 99.46401
## 1738 47362 11.6450321 sc 34.29054
## 1739 47375 0.7514160 sc 98.22866
## 1740 47381 0.7072351 sc 94.78186
## 1741 47382 1.5671769 sc 99.62035library(patchwork)
p1 <- ggplot(df_long,aes(Dnsty,y,col=凡例))+
geom_point(colour="blue")+
labs(title = "散布図(1):",subtitle = "可住地人口密度と自地域就業率(全国)",x="可住地人口密度",y="自地域就業率",caption = "2015年")+
theme_minimal(base_family="HiraKakuProN-W3")p2 <- ggplot(df_long,aes(Dnsty,y,col=凡例))+
geom_point(colour="seagreen")+
labs(title = "散布図(2)",subtitle = "x: 可住地人口密度\ny: 自地域就業率",x="x",caption = "2015年")+
theme_minimal(base_family="HiraKakuProN-W3")p3 <- ggplot(df_long,aes(Dnsty,y,col=凡例))+
geom_point(colour="skyblue")+
labs(title = "散布図(3)",subtitle = paste("x=",colnames(dat)[26],"\ny=",colnames(dat)[30]),x="x",caption = "2015年")+
theme_minimal(base_family="HiraKakuProN-W3")
散布図(1)は、できるだけわかりやすく説明を日本語で入れるもの。散布図(3)は、できるだけ自動的に図を作成する場合の雛形。サブタイトルの部分が自動的に挿入される。データフレームの変数名が反映する。
次の散布図(4)のコードと実行結果はつぎの通り。
df_long <- gather(dat[,c(2,31,33)],"凡例","y",-Code,-r_nuclear_relativeH)
ggplot(df_long,aes(r_nuclear_relativeH,y,col=凡例))+
geom_point(colour="skyblue")+
ylim(0,20)+
labs(title = "散布図(4)",subtitle = paste("x=",colnames(dat)[33],"\ny=",colnames(dat)[31]),x="x",caption = "2015年")+
theme_minimal(base_family="HiraKakuProN-W3")
下のボールドの部分を書き換えると同じようなグラフが描ける。取りあげる2変数の列番号、x軸に指定する変数の名前、サブタイトル用の変数番号の部分である。
df_long <- gather(dat[,c(2,31,33)],“凡例”,“y”,-Code,-r_nuclear_relativeH)
ggplot(df_long,aes(r_nuclear_relativeH,y,col=凡例))+
geom_point(colour=“skyblue”)+
ylim(0,20)+ labs(title = “散布図(4)”,subtitle =
paste(“x=”,colnames(dat)[33],“=”,colnames(dat)[31]),x=“x”,caption = “2015年”)+
theme_minimal(base_family=“HiraKakuProN-W3”)
2.4 福岡県
## Warning: ggrepel: 20 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: ggrepel: 34 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: ggrepel: 25 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: ggrepel: 14 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: ggrepel: 13 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
2.5 群馬県
## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
3. 主成分分析
福岡県内各市町村単位のデータを用いて分析をおこなう。 データは、RDataファイルではなく、RDSファイル2として保存されている。
3.1 準備
RDSファイルから読み込んだデータの全項目のリストはつぎの通りである。 各列の1行目に項目の解説の行が含まれている。（その行は後で削除する。）
## X1 AREA Code
## X2 AREA
## A1101 A1101_Total population (Both sexes)[person]
## A5101 A5101_Number of in-migrants from other municipalities[person]
## A5102 A5102_Number of out-migrants to other municipalities[person]
## A6108 A6108_Rate of day to night population[%]
## B1103 B1103_Inhabitable area[ha]
## D2201 D2201_Financial strength index (municipalities)[-]
## F1102 F1102_(Population census) Employed persons[person]
## F2701 F2701_Population working in the same municipalities[person]
## A7101 A7101_Number of households (Total)[households]
## A710101 A710101_Number of private households[households]
## A710102 A710102_Number of institutional households[households]
## A710201 A710201_No. of private household members[person]
## A710202 A710202_No. of institutional household members[person]
## A7103 A7103_No. of basic resident register households (Japanese)[households]
## A810101 A810101_No. of relatives households[households]
## A810102 A810102_Nuclear families[households]
## A810103 A810103_Relatives households excluding nuclear families[households]
## A810104 A810104_Households including non-relatives[households]
## A810105 A810105_One-person households[households]
## A811102 A811102_No. of nuclear families with household members 65 years of age and over[households]
## A8201 A8201_No. of aged-couple households (Only aged couple)[households]
## A8301 A8301_No. of aged-single-person households (Aged 65 and over)[households]
## A830101 A830101_No. of aged-single-person households (Aged 65 and over, male)[households]
## A830102 A830102_No. of aged-single-person households (Aged 65 and over, female)[households]
## A8401 A8401_No. of mother-child(ren) households[households]
## A840101 A840101_No. of mother-child households (1 child)[households]
## A840102 A840102_No. of mother-children households (2 children)[households]
## A840103 A840103_No. of mother-children households (3 children or more)[households]
## A8501 A8501_No. of father-child(ren) households[households]
## A850101 A850101_No. of father-child households (1 child)[households]
## A850102 A850102_No. of father-children households (2 children)[households]
## A850103 A850103_No. of father-children households (3 children or more)[households]
## A110101 A110101_Total population (Male)[person]
## A110102 A110102_Total population (Female)[person]
## A1301 A1301_Total population (0-14)[person]
## A1302 A1302_Total population (15-64)[person]
## A1303 A1303_Total population (65 and over)[person]
## A1419 A1419_Population (75 and over, total)[person]
## A141901 A141901_Population (75 and over, male)[person]
## A141902 A141902_Population (75 and over, female)[person]
## A161001 A161001_Ratio of never married persons (population 15 years old and over)[%]
## A1700 A1700_Foreigners (Total)[person]
## A4101 A4101_Live births[person]
## A4200 A4200_Deaths[person]
## A9101 A9101_Marriages[pairs]
## A9201 A9201_Divorces[pairs]
## D2205 D2205_Basic financial revenues (local public finance, prefecture)[thousand yen]
## D2206 D2206_Basic financial needs (local public finance, prefectures)[thousand yen]3.1.1 datの作成
特定の分析に利用する列を取り出す。データ名をdatとする。
式で新変数を構成する。 列を追加するためにmutate()を使う。 基本的に変数の順序は変更しない。
列番号と変数名の対応関係
# 列番号と変数名の対応関係の確認（オブジェクト名は「dat」）
for (i in 1:ncol(dat)) {
cat(i,"===",colnames(dat)[i],"\n")
}## 1 === X1
## 2 === X2
## 3 === A1101
## 4 === A6108
## 5 === B1103
## 6 === A5101
## 7 === A5102
## 8 === D2201
## 9 === F1102
## 10 === F2701
## 11 === A710101
## 12 === A710201
## 13 === A710202
## 14 === A810101
## 15 === A810102
## 16 === A810103
## 17 === A810104
## 18 === A810105
## 19 === A1303
## 20 === A4101
## 21 === A4200
## 22 === D2205
## 23 === D2206
## 24 === A7101
## 25 === Dnsty
## 26 === In
## 27 === Out
## 28 === OutOverIn
## 29 === sc
## 30 === r_nuclear
## 31 === r_nuclear_relativeH
## 32 === r_OnePerson
## 33 === r_instHm
## 34 === aged
## 35 === naturalIncreae相関行列
# 相関行列（全変数）
round(cor(dat[,c(3:ncol(dat))]),digit=2)## A1101 A6108 B1103 A5101 A5102 D2201 F1102 F2701 A710101
## A1101 1.00 0.25 0.76 0.97 0.98 0.35 1.00 0.99 1.00
## A6108 0.25 1.00 0.32 0.23 0.24 0.36 0.25 0.27 0.25
## B1103 0.76 0.32 1.00 0.64 0.67 0.23 0.77 0.83 0.73
## A5101 0.97 0.23 0.64 1.00 1.00 0.34 0.97 0.94 0.99
## A5102 0.98 0.24 0.67 1.00 1.00 0.35 0.98 0.95 0.99
## D2201 0.35 0.36 0.23 0.34 0.35 1.00 0.36 0.34 0.34
## F1102 1.00 0.25 0.77 0.97 0.98 0.36 1.00 0.99 1.00
## F2701 0.99 0.27 0.83 0.94 0.95 0.34 0.99 1.00 0.98
## A710101 1.00 0.25 0.73 0.99 0.99 0.34 1.00 0.98 1.00
## A710201 1.00 0.25 0.76 0.97 0.99 0.35 1.00 0.99 1.00
## A710202 0.98 0.27 0.82 0.92 0.94 0.32 0.98 0.99 0.97
## A810101 1.00 0.24 0.78 0.96 0.97 0.36 1.00 0.99 0.99
## A810102 1.00 0.24 0.76 0.96 0.98 0.36 1.00 0.99 0.99
## A810103 0.96 0.27 0.89 0.88 0.90 0.34 0.96 0.99 0.94
## A810104 0.99 0.24 0.69 0.99 1.00 0.34 0.99 0.96 1.00
## A810105 0.98 0.25 0.67 1.00 1.00 0.32 0.98 0.95 0.99
## A1303 0.98 0.25 0.82 0.92 0.94 0.33 0.98 0.99 0.97
## A4101 1.00 0.24 0.73 0.99 0.99 0.37 1.00 0.98 1.00
## A4200 0.98 0.27 0.83 0.90 0.92 0.32 0.98 0.99 0.96
## D2205 1.00 0.26 0.73 0.98 0.99 0.35 1.00 0.98 1.00
## D2206 1.00 0.26 0.79 0.96 0.97 0.33 1.00 0.99 0.99
## A7101 1.00 0.25 0.73 0.99 0.99 0.34 1.00 0.98 1.00
## Dnsty 0.47 -0.13 0.10 0.49 0.49 0.55 0.47 0.41 0.47
## In 0.28 0.00 -0.09 0.33 0.32 0.63 0.28 0.22 0.28
## Out 0.27 -0.07 -0.12 0.33 0.32 0.43 0.27 0.21 0.28
## OutOverIn -0.14 -0.04 0.02 -0.16 -0.15 -0.57 -0.15 -0.12 -0.14
## sc 0.16 0.51 0.57 0.08 0.10 -0.09 0.16 0.24 0.14
## r_nuclear -0.36 -0.42 -0.42 -0.35 -0.36 0.30 -0.36 -0.39 -0.37
## r_nuclear_relativeH 0.25 -0.08 -0.04 0.25 0.25 0.47 0.25 0.22 0.25
## r_OnePerson 0.56 0.35 0.36 0.55 0.55 0.09 0.55 0.55 0.56
## r_instHm -0.18 -0.01 -0.12 -0.18 -0.18 -0.66 -0.19 -0.17 -0.17
## aged -0.25 0.04 -0.06 -0.27 -0.27 -0.82 -0.26 -0.22 -0.25
## naturalIncreae 0.37 0.00 -0.09 0.56 0.52 0.24 0.37 0.26 0.42
## A710201 A710202 A810101 A810102 A810103 A810104 A810105
## A1101 1.00 0.98 1.00 1.00 0.96 0.99 0.98
## A6108 0.25 0.27 0.24 0.24 0.27 0.24 0.25
## B1103 0.76 0.82 0.78 0.76 0.89 0.69 0.67
## A5101 0.97 0.92 0.96 0.96 0.88 0.99 1.00
## A5102 0.99 0.94 0.97 0.98 0.90 1.00 1.00
## D2201 0.35 0.32 0.36 0.36 0.34 0.34 0.32
## F1102 1.00 0.98 1.00 1.00 0.96 0.99 0.98
## F2701 0.99 0.99 0.99 0.99 0.99 0.96 0.95
## A710101 1.00 0.97 0.99 0.99 0.94 1.00 0.99
## A710201 1.00 0.98 1.00 1.00 0.96 0.99 0.98
## A710202 0.98 1.00 0.99 0.99 0.98 0.95 0.94
## A810101 1.00 0.99 1.00 1.00 0.97 0.98 0.97
## A810102 1.00 0.99 1.00 1.00 0.96 0.98 0.97
## A810103 0.96 0.98 0.97 0.96 1.00 0.92 0.90
## A810104 0.99 0.95 0.98 0.98 0.92 1.00 1.00
## A810105 0.98 0.94 0.97 0.97 0.90 1.00 1.00
## A1303 0.98 1.00 0.99 0.99 0.98 0.95 0.94
## A4101 1.00 0.97 0.99 0.99 0.94 1.00 0.99
## A4200 0.98 1.00 0.99 0.98 0.98 0.94 0.92
## D2205 1.00 0.98 0.99 1.00 0.94 0.99 0.99
## D2206 1.00 0.99 1.00 1.00 0.97 0.98 0.97
## A7101 1.00 0.97 0.99 1.00 0.94 1.00 0.99
## Dnsty 0.47 0.41 0.47 0.47 0.38 0.48 0.46
## In 0.28 0.21 0.27 0.28 0.15 0.31 0.29
## Out 0.27 0.21 0.26 0.28 0.14 0.31 0.29
## OutOverIn -0.15 -0.12 -0.15 -0.15 -0.10 -0.15 -0.13
## sc 0.16 0.23 0.17 0.15 0.32 0.11 0.12
## r_nuclear -0.36 -0.37 -0.35 -0.34 -0.39 -0.36 -0.39
## r_nuclear_relativeH 0.25 0.24 0.26 0.27 0.15 0.26 0.24
## r_OnePerson 0.56 0.56 0.55 0.55 0.50 0.56 0.57
## r_instHm -0.18 -0.12 -0.19 -0.18 -0.20 -0.18 -0.15
## aged -0.26 -0.20 -0.26 -0.26 -0.22 -0.27 -0.23
## naturalIncreae 0.38 0.21 0.33 0.34 0.17 0.49 0.51
## A1303 A4101 A4200 D2205 D2206 A7101 Dnsty In Out
## A1101 0.98 1.00 0.98 1.00 1.00 1.00 0.47 0.28 0.27
## A6108 0.25 0.24 0.27 0.26 0.26 0.25 -0.13 0.00 -0.07
## B1103 0.82 0.73 0.83 0.73 0.79 0.73 0.10 -0.09 -0.12
## A5101 0.92 0.99 0.90 0.98 0.96 0.99 0.49 0.33 0.33
## A5102 0.94 0.99 0.92 0.99 0.97 0.99 0.49 0.32 0.32
## D2201 0.33 0.37 0.32 0.35 0.33 0.34 0.55 0.63 0.43
## F1102 0.98 1.00 0.98 1.00 1.00 1.00 0.47 0.28 0.27
## F2701 0.99 0.98 0.99 0.98 0.99 0.98 0.41 0.22 0.21
## A710101 0.97 1.00 0.96 1.00 0.99 1.00 0.47 0.28 0.28
## A710201 0.98 1.00 0.98 1.00 1.00 1.00 0.47 0.28 0.27
## A710202 1.00 0.97 1.00 0.98 0.99 0.97 0.41 0.21 0.21
## A810101 0.99 0.99 0.99 0.99 1.00 0.99 0.47 0.27 0.26
## A810102 0.99 0.99 0.98 1.00 1.00 1.00 0.47 0.28 0.28
## A810103 0.98 0.94 0.98 0.94 0.97 0.94 0.38 0.15 0.14
## A810104 0.95 1.00 0.94 0.99 0.98 1.00 0.48 0.31 0.31
## A810105 0.94 0.99 0.92 0.99 0.97 0.99 0.46 0.29 0.29
## A1303 1.00 0.97 1.00 0.98 0.99 0.97 0.42 0.22 0.21
## A4101 0.97 1.00 0.96 1.00 0.99 1.00 0.49 0.31 0.30
## A4200 1.00 0.96 1.00 0.97 0.99 0.96 0.40 0.19 0.19
## D2205 0.98 1.00 0.97 1.00 0.99 1.00 0.46 0.28 0.28
## D2206 0.99 0.99 0.99 0.99 1.00 0.99 0.43 0.24 0.24
## A7101 0.97 1.00 0.96 1.00 0.99 1.00 0.47 0.28 0.28
## Dnsty 0.42 0.49 0.40 0.46 0.43 0.47 1.00 0.62 0.65
## In 0.22 0.31 0.19 0.28 0.24 0.28 0.62 1.00 0.81
## Out 0.21 0.30 0.19 0.28 0.24 0.28 0.65 0.81 1.00
## OutOverIn -0.12 -0.16 -0.10 -0.14 -0.12 -0.14 -0.30 -0.70 -0.22
## sc 0.21 0.14 0.24 0.14 0.19 0.14 -0.34 -0.50 -0.46
## r_nuclear -0.36 -0.35 -0.37 -0.37 -0.38 -0.37 0.21 0.41 0.22
## r_nuclear_relativeH 0.24 0.26 0.23 0.24 0.23 0.25 0.58 0.64 0.60
## r_OnePerson 0.54 0.56 0.55 0.55 0.56 0.56 0.28 0.14 0.30
## r_instHm -0.17 -0.19 -0.16 -0.17 -0.16 -0.17 -0.42 -0.36 -0.20
## aged -0.21 -0.28 -0.19 -0.24 -0.22 -0.25 -0.65 -0.76 -0.61
## naturalIncreae 0.20 0.44 0.16 0.40 0.31 0.42 0.41 0.46 0.43
## OutOverIn sc r_nuclear r_nuclear_relativeH r_OnePerson
## A1101 -0.14 0.16 -0.36 0.25 0.56
## A6108 -0.04 0.51 -0.42 -0.08 0.35
## B1103 0.02 0.57 -0.42 -0.04 0.36
## A5101 -0.16 0.08 -0.35 0.25 0.55
## A5102 -0.15 0.10 -0.36 0.25 0.55
## D2201 -0.57 -0.09 0.30 0.47 0.09
## F1102 -0.15 0.16 -0.36 0.25 0.55
## F2701 -0.12 0.24 -0.39 0.22 0.55
## A710101 -0.14 0.14 -0.37 0.25 0.56
## A710201 -0.15 0.16 -0.36 0.25 0.56
## A710202 -0.12 0.23 -0.37 0.24 0.56
## A810101 -0.15 0.17 -0.35 0.26 0.55
## A810102 -0.15 0.15 -0.34 0.27 0.55
## A810103 -0.10 0.32 -0.39 0.15 0.50
## A810104 -0.15 0.11 -0.36 0.26 0.56
## A810105 -0.13 0.12 -0.39 0.24 0.57
## A1303 -0.12 0.21 -0.36 0.24 0.54
## A4101 -0.16 0.14 -0.35 0.26 0.56
## A4200 -0.10 0.24 -0.37 0.23 0.55
## D2205 -0.14 0.14 -0.37 0.24 0.55
## D2206 -0.12 0.19 -0.38 0.23 0.56
## A7101 -0.14 0.14 -0.37 0.25 0.56
## Dnsty -0.30 -0.34 0.21 0.58 0.28
## In -0.70 -0.50 0.41 0.64 0.14
## Out -0.22 -0.46 0.22 0.60 0.30
## OutOverIn 1.00 0.36 -0.49 -0.44 0.11
## sc 0.36 1.00 -0.62 -0.47 0.20
## r_nuclear -0.49 -0.62 1.00 0.46 -0.59
## r_nuclear_relativeH -0.44 -0.47 0.46 1.00 0.44
## r_OnePerson 0.11 0.20 -0.59 0.44 1.00
## r_instHm 0.39 0.16 -0.26 -0.18 0.12
## aged 0.56 0.32 -0.47 -0.48 0.08
## naturalIncreae -0.23 -0.29 -0.06 0.20 0.21
## r_instHm aged naturalIncreae
## A1101 -0.18 -0.25 0.37
## A6108 -0.01 0.04 0.00
## B1103 -0.12 -0.06 -0.09
## A5101 -0.18 -0.27 0.56
## A5102 -0.18 -0.27 0.52
## D2201 -0.66 -0.82 0.24
## F1102 -0.19 -0.26 0.37
## F2701 -0.17 -0.22 0.26
## A710101 -0.17 -0.25 0.42
## A710201 -0.18 -0.26 0.38
## A710202 -0.12 -0.20 0.21
## A810101 -0.19 -0.26 0.33
## A810102 -0.18 -0.26 0.34
## A810103 -0.20 -0.22 0.17
## A810104 -0.18 -0.27 0.49
## A810105 -0.15 -0.23 0.51
## A1303 -0.17 -0.21 0.20
## A4101 -0.19 -0.28 0.44
## A4200 -0.16 -0.19 0.16
## D2205 -0.17 -0.24 0.40
## D2206 -0.16 -0.22 0.31
## A7101 -0.17 -0.25 0.42
## Dnsty -0.42 -0.65 0.41
## In -0.36 -0.76 0.46
## Out -0.20 -0.61 0.43
## OutOverIn 0.39 0.56 -0.23
## sc 0.16 0.32 -0.29
## r_nuclear -0.26 -0.47 -0.06
## r_nuclear_relativeH -0.18 -0.48 0.20
## r_OnePerson 0.12 0.08 0.21
## r_instHm 1.00 0.59 -0.17
## aged 0.59 1.00 -0.36
## naturalIncreae -0.17 -0.36 1.00# 相関図行列（例）
pairs(dat[,c(5,25,30,32,33,34,35,26,28,29,4)])
3.1.2 dat_sの作成
以下のような変数が含まれているデータファイルを分析の対象とする。 元のデータdatから、この12変数のみを取り出し、dat_sと名付ける。
## 1 === X2 地域名
## 2 === B11 可住地面積
## 3 === Dns 可住地人口密度
## 4 === r_n 核家族世帯率
## 5 === r_O 単独世帯率
## 6 === r_H 施設居住者率
## 7 === agd 高齢者人口比率
## 8 === ntI 人口の自然増加
## 9 === In 転入者率
## 10 === OOI 転出超過率
## 11 === sc 自地域就業率
## 12 === A61 昼夜間人口比率3.2 princompを使う
まず、データを確認してみる。1列目は地名の略号である。
nrow(dat_s)## [1] 60row.names(dat_s) <- 1:60dat_s[1:10,] # データの確認## X2 B11 Dns r_n r_O r_H agd ntI In
## 1 Ktky 29609 32.466007 56.09032 0.3700863 2.921295 28.82805 -2640 0.04197606
## 2 Fkks 23174 66.396867 45.68579 0.4968409 2.138065 20.29862 3631 0.07187390
## 3 Omts 6392 18.360451 54.89949 0.3417056 4.049932 34.54499 -965 0.02481254
## 4 Krms 19530 15.594060 55.10617 0.3289765 2.824805 25.03875 -148 0.03757322
## 5 Ngts 3992 14.315130 59.97760 0.3029677 3.627550 31.39502 -201 0.03559304
## 6 Izks 10746 12.018053 54.80681 0.3580156 3.909529 28.81235 -434 0.03371378
## 7 Tgws 3900 12.420769 54.98659 0.3743054 5.090729 31.95640 -302 0.03707603
## 8 Yngw 7715 8.785094 55.39362 0.2330541 2.878558 30.71101 -386 0.02236747
## 9 Ymsh 16535 3.895252 55.26089 0.2211979 3.970004 33.30487 -503 0.02280773
## 10 Chkg 4171 11.589307 61.03230 0.2362745 2.426612 25.76801 -53 0.03963673
## OOI sc A61
## 1 1.0765285 57.24755 102.3
## 2 0.9197041 46.12117 110.8
## 3 1.1895604 73.55213 104.5
## 4 0.9614612 67.74829 99.5
## 5 1.0294985 51.88872 105.1
## 6 1.0569591 67.93815 101.5
## 7 1.0579065 59.17267 109.0
## 8 1.3449868 57.89591 92.4
## 9 1.3689585 71.38323 102.6
## 10 0.9968685 45.77679 96.4# 主成分分析(princomp)
# dat_sの1列目の地名を削除して計算
d <- dat_s[,-1]
model <- princomp(d,cor = TRUE)
# 結果の出力
summary(model,loading=TRUE)## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## Standard deviation 2.0122464 1.5548599 1.1441406 0.94788027 0.82580397
## Proportion of Variance 0.3681032 0.2197809 0.1190053 0.08167973 0.06199565
## Cumulative Proportion 0.3681032 0.5878841 0.7068893 0.78856907 0.85056472
## Comp.6 Comp.7 Comp.8 Comp.9 Comp.10
## Standard deviation 0.72154851 0.65022834 0.47948327 0.46862734 0.38673934
## Proportion of Variance 0.04733021 0.03843608 0.02090038 0.01996469 0.01359703
## Cumulative Proportion 0.89789493 0.93633101 0.95723139 0.97719608 0.99079311
## Comp.11
## Standard deviation 0.318238561
## Proportion of Variance 0.009206889
## Cumulative Proportion 1.000000000
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9 Comp.10
## B11 0.131 0.436 0.335 0.304 0.261 0.510 0.434 0.232
## Dns -0.332 0.271 -0.182 0.421 0.167 -0.209 -0.289 0.171 -0.621 0.189
## r_n -0.345 -0.354 0.230 0.194 -0.223 0.330 -0.664
## r_O 0.472 -0.406 0.359 -0.275 0.293 -0.107 0.143 -0.522
## r_H 0.279 -0.163 -0.390 -0.312 0.480 0.364 -0.493
## agd 0.412 -0.189 -0.198 -0.102 0.426 0.280 -0.398
## ntI -0.222 0.224 -0.468 -0.113 -0.606 0.469 0.154 -0.220
## In -0.421 0.194 -0.243 0.222 -0.277 0.284 0.292
## OOI 0.356 -0.246 0.449 -0.213 -0.300 -0.432 0.133 0.427
## sc 0.351 0.275 0.335 -0.116 0.157 -0.297 -0.536 -0.335 -0.308
## A61 0.144 0.394 0.214 -0.588 -0.168 -0.392 -0.255 0.421
## Comp.11
## B11
## Dns
## r_n 0.250
## r_O
## r_H -0.175
## agd 0.565
## ntI
## In 0.647
## OOI 0.286
## sc 0.258
## A61# 寄与率
#((model$sdev[1])^2)/ncol(d)
#((model$sdev[2])^2)/ncol(d)
#((model$sdev[3])^2)/ncol(d)各個体の主成分得点
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6
## 1 1.880 3.092 3.124 1.892 3.855 -0.590
## 2 -2.594 8.068 -3.322 0.317 -1.525 2.206
## 3 2.950 1.312 0.705 0.129 0.572 -0.922
## 4 0.973 2.923 1.639 0.501 0.421 1.189
## 5 0.737 0.611 0.309 -0.824 0.019 -0.465
## 6 1.832 2.046 0.677 -0.109 0.590 0.115
## 7 1.827 1.561 -0.481 -1.272 0.682 -0.581
## 8 2.057 -0.052 0.809 1.186 -0.979 -0.054
## 9 3.343 0.876 1.704 0.830 -0.502 0.824
## 10 -0.544 -0.010 0.940 -0.282 -0.654 -0.002
## 11 1.764 0.474 0.532 0.104 -1.216 -0.729
## 12 -0.156 0.118 0.806 0.125 -0.140 0.093
## 13 2.391 -0.184 -0.041 -0.705 0.003 0.068
## 14 -0.313 -0.650 -0.498 0.634 0.142 -0.853
## 15 -1.292 -1.249 0.537 0.354 0.064 0.358
## 16 -1.937 -0.182 0.271 0.185 0.260 0.085
## 17 -4.525 0.717 -0.752 2.025 0.566 -0.794
## 18 -3.527 0.723 -0.908 1.237 0.074 -0.636
## 19 -0.713 0.073 0.383 0.652 0.042 0.395
## 20 -2.387 0.291 -0.589 0.016 0.477 -0.659
## 21 -1.413 0.039 0.533 -0.191 -0.289 -0.311
## 22 -2.420 -1.070 0.846 -1.115 0.837 1.147
## 23 2.030 -0.220 0.856 0.991 -1.506 0.085
## 24 2.338 1.158 0.197 -2.110 -0.006 -0.602
## 25 2.466 -0.325 -0.729 0.422 0.387 -0.098
## 26 2.593 1.509 1.491 0.432 -0.897 -0.101
## 27 1.868 -0.822 0.741 0.790 -0.724 0.739
## 28 -0.574 -0.352 1.896 0.751 -0.013 1.346
## 29 -2.850 -0.895 0.959 0.809 -0.713 -0.164
## 30 -1.192 -1.558 -0.362 0.231 0.562 0.500
## 31 -2.289 -1.089 0.190 0.302 0.045 0.227
## 32 -2.482 0.046 -0.471 1.625 -0.237 -1.129
## 33 -2.215 -0.502 0.154 -0.431 0.204 0.085
## 34 -5.107 0.258 1.369 -2.311 0.567 0.280
## 35 -0.623 0.430 1.587 -2.828 -1.272 -1.197
## 36 -3.759 1.128 0.053 0.115 -0.403 -0.774
## 37 -0.418 -0.283 -1.423 -0.734 1.201 0.120
## 38 -1.164 -0.455 -0.846 0.844 0.009 -0.911
## 39 -0.877 -2.024 0.097 0.120 0.390 0.680
## 40 -1.140 -1.380 0.708 -0.293 -0.534 -0.259
## 41 0.744 -0.016 -1.550 -1.340 0.626 -0.410
## 42 0.717 -0.568 -0.003 -0.018 -0.664 -0.784
## 43 -0.798 -1.195 -0.557 0.112 0.081 0.188
## 44 -0.660 -1.675 0.470 -0.142 0.053 1.312
## 45 3.325 -0.414 -0.567 -0.028 -1.210 0.470
## 46 -0.273 -1.417 0.008 -0.041 -0.406 0.633
## 47 -0.793 -1.756 0.782 0.764 -1.258 0.204
## 48 -0.303 0.104 0.108 -1.343 0.063 0.553
## 49 1.021 -0.699 -1.163 -0.500 0.219 -0.249
## 50 2.143 -0.912 -1.323 0.211 0.132 0.179
## 51 0.156 -0.674 -2.207 -0.008 1.058 -0.042
## 52 2.630 0.173 -2.042 1.324 -0.451 -1.353
## 53 1.368 -0.804 -2.143 -0.894 1.071 0.488
## 54 2.503 -1.703 -1.755 0.599 -0.234 0.065
## 55 1.926 -1.056 -1.823 -0.008 0.916 0.377
## 56 -0.151 2.981 0.514 -1.813 -0.444 -1.519
## 57 1.433 -1.484 -0.114 -0.582 0.436 0.846
## 58 -0.366 -0.959 -0.360 0.346 -0.646 -0.801
## 59 -0.446 -1.422 0.664 -0.456 -0.504 0.344
## 60 1.289 -0.654 -0.630 -0.601 0.802 0.7883.2.1 個体のプロット
3.2.2 変数のプロット
princompの計算結果がmodelというオブジェクト（データフレーム）に保存されている。
これから数値を取り出して各変数の「負荷量(loadings)」を2次元の図に描く。（この負荷量の図は、主成分負荷量の図と形状は同じであるが、同じではない。後者は、主成分と各変数の間の相関係数である。）
最初の行mdl_ldgs <- model$loadings[,c(1,2)]
の最後の数字を入れ替えることで第1主成分と第2主成分、第1主成分と第3主成分のグラフを作成できる。(ラベルは変更する必要がある。)
ggplot2には、long formatのデータを渡している。
mdl_ldgs <- model$loadings[,c(1,2)]
mdl_ldgs <- as.data.frame(mdl_ldgs)
mdl_long <- gather(mdl_ldgs,"key","value",-Comp.1)
rn <- rownames(mdl_ldgs)
l1 <- round(((model$sdev[1])^2)/ncol(d)*100,digits = 2)
l2 <- round(((model$sdev[2])^2)/ncol(d)*100,digits = 2)
l3 <- round(((model$sdev[3])^2)/ncol(d)*100,digits = 2)library(ggrepel)
library(stringr)
library(patchwork)
p1 <- ggplot(mdl_long,aes(Comp.1,value,col=key))+
geom_point(col="blue")+
labs(x=str_c("PC1(",l1,"%)"),y=str_c("PC2","(",l2,"%)"),title = "Fukuoka")+
geom_text_repel(aes(label=rn),size=3) +
geom_hline(yintercept = 0,linetype="dashed")+
geom_vline(xintercept = 0,linetype="dashed")+
annotate("text",x=0.30,y=0.20,label="sc(0.35, 0.28)")+
theme_minimal(base_family="HiraKakuProN-W3")+
theme(legend.position="none")p1 + p2
3.3 描画にggplot2を使う
個体の主成分スコアのプロット。
m1 <- model$scores[,1]
m2 <- model$scores[,2]
m1m2 <- data.frame(pc1 =m1,pc2=m2)## pc1 pc2
## Min. :-5.1066 Min. :-2.0244
## 1st Qu.:-1.1713 1st Qu.:-0.9236
## Median :-0.2878 Median :-0.2517
## Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 1.8412 3rd Qu.: 0.4408
## Max. : 3.3426 Max. : 8.0680図(-6.1, 6.1)
図(-2.5, 2.5)
別の方法。 long formatを以下のように利用。
mdl_scrs <- model$scores[,c(1,2)]
mdl_scrs <- as.data.frame(mdl_scrs)
mdl_long <- gather(mdl_scrs,"key","value",-Comp.1)mdl_scrs <- model$scores[,c(1,2)]
mdl_scrs <- as.data.frame(mdl_scrs)
mdl_long <- gather(mdl_scrs,"key","value",-Comp.1)
rn <- rownames(mdl_scrs) # 地域番号で表記
#rn <- abbreviate(dat_s$X2)
library(ggrepel)
ggplot(mdl_long,aes(Comp.1,value,col=key))+
geom_point(col="black")+ # 点の色は変更可能
labs(x=str_c("PC1(",l1,"%)"),y=str_c("PC2","(",l2,"%)"),title = "Fukuoka")+
geom_text_repel(aes(label=rn),size=3,col="black") + # 番号の色は変更可能
geom_hline(yintercept = 0,linetype="dashed")+
geom_vline(xintercept = 0,linetype="dashed")+
theme_minimal(base_family="HiraKakuProN-W3")+
theme(legend.position="none")
3.4 FactoMineRを使う
主成分負荷量が簡単に取り出せるので便利である。
## eigenvalue percentage of variance cumulative percentage of variance
## comp 1 4.0491355 36.8103227 36.81032
## comp 2 2.4175895 21.9780860 58.78841
## comp 3 1.3090578 11.9005254 70.68893
## comp 4 0.8984770 8.1679727 78.85691
## comp 5 0.6819522 6.1995654 85.05647
## comp 6 0.5206323 4.7330205 89.78949
## comp 7 0.4227969 3.8436081 93.63310
## comp 8 0.2299042 2.0900382 95.72314
## comp 9 0.2196116 1.9964690 97.71961
## comp 10 0.1495673 1.3597029 99.07931
## comp 11 0.1012758 0.9206889 100.00000## Dim.1 Dim.2 Dim.3
## B11 0.2630344 0.6785060 -0.3828134
## Dns -0.6677016 0.4211047 0.2085956
## r_n -0.6945110 -0.5496752 -0.2625984
## r_O 0.1819990 0.7336440 0.4648059
## r_H 0.5619070 -0.2529011 0.4463726
## agd 0.8297779 -0.2940360 0.2263402
## ntI -0.4476360 0.3485525 0.5356338
## In -0.8475153 0.3022256 0.1013067
## OOI 0.7163329 -0.1315104 0.2816353
## sc 0.7059704 0.4277115 -0.3833217
## A61 0.2901654 0.6132178 -0.2451493
# 第1から第3主成分までの主成分負荷量の散布図
pairs(dat.PCA$var$cor[,1:3],xlim=c(-1,1),ylim=c(-1,1))
以下の図では、図中に変数scの座標を入れてみた。また、寄与率も各軸に表示した。
par(pty="s",family="HiraKakuProN-W3",xpd=TRUE)
plot(dat.PCA$var$cor[,1],dat.PCA$var$cor[,2],type="n",main="各変数の主成分負荷量",xlab = str_c("PC1(",l1,"%)"), ylab = str_c("PC2(",l2,"%)") ,cex.main=0.8,cex.axis=0.8,cex.lab=0.8)
text(dat.PCA$var$cor[,1],dat.PCA$var$cor[,2],rownames(dat.PCA$var$coord),cex=0.8)
#abline(h=0,lty=3)
#abline(v=0,lty=3)
text(x=0.65,y=0.34,"(0.71,-0.13)",cex=0.7,col="blue")
第1主成分と第3主成分のグラフ
第2主成分と第3主成分のグラフ
3.5 prcompを使う
prcompの説明：
## Default S3 method:
prcomp(x, retx = TRUE, center = TRUE, scale. = FALSE,
tol = NULL, rank. = NULL, ...)分析対象オブジェクトをdとし、主成分分析をおこない、結果をdat.prcomに保存する。
scale=TRUEとすること。
# 分析対象オブジェクトをdとする。
d <- dat_s[,-1] # 1列目（地域名）を外す
dat.prcom <- prcomp(d,scale = TRUE)
summary(dat.prcom) # Importance of componetsの表## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.0122 1.5549 1.1441 0.94788 0.8258 0.72155 0.65023
## Proportion of Variance 0.3681 0.2198 0.1190 0.08168 0.0620 0.04733 0.03844
## Cumulative Proportion 0.3681 0.5879 0.7069 0.78857 0.8506 0.89789 0.93633
## PC8 PC9 PC10 PC11
## Standard deviation 0.4795 0.46863 0.3867 0.31824
## Proportion of Variance 0.0209 0.01996 0.0136 0.00921
## Cumulative Proportion 0.9572 0.97720 0.9908 1.00000#dat.prcom # Rotationの表Rotationの表に出てくる数値は、エヴェリットの翻訳に出てくる「負荷量」と同じ。
主成分負荷量を計算するために、 関数prcom33を使う。
prcomp3 <- function(df) {
df<- na.omit(df)
ans <- prcomp(df,scale=TRUE)
ans$loadings <- t(t(ans$rotation)*ans$sdev)
ans$eigenvalues <- ans$sdev^2
print.default(ans)
invisible(ans)
# Aoki
}上記関数を使用し、 結果をオブジェクトans3に「付値」する。 最後の方のloadingsの表が、主成分負荷量の数値。
biplot(ans3,choices = c(1,2))
biplot(ans3,choices = c(1,3))
biplot。矢印の先が「主成分負荷量」の値である。各ケースの「主成分得点」は番号の位置で表される。
plot(ans3$loadings[,1],ans3$loadings[,2],asp=1)
abline(h=0,v=0)
text(ans3$loadings[,1],ans3$loadings[,2],
labels=row.names(ans3$loadings),pos=3)
library(tidyverse)
mdl_ldgs <- ans3$loadings[,c(1,2)]
mdl_ldgs <- as.data.frame(mdl_ldgs)
mdl_long <- gather(mdl_ldgs,"key","value",-PC1)
rn <- rownames(mdl_ldgs)
l1 <- round(((ans3$sdev[1])^2)/ncol(d)*100,digits = 2)
l2 <- round(((ans3$sdev[2])^2)/ncol(d)*100,digits = 2)
l3 <- round(((ans3$sdev[3])^2)/ncol(d)*100,digits = 2)
library(stringr)
library(ggrepel)
ggplot(mdl_long,aes(PC1,value,col=key))+
geom_point(col="blue")+
labs(x=str_c("PC1(",l1,"%)"),y=str_c("PC2","(",l2,"%)"),title = "Fukuoka")+
geom_text_repel(aes(label=rn),size=3,col="Black") +
geom_hline(yintercept = 0,linetype="dashed")+
geom_vline(xintercept = 0,linetype="dashed")+
xlim(-1,1)+
ylim(-1,1)+
theme_minimal(base_family="HiraKakuProN-W3")+
theme(legend.position="none")
4. 対応分析を行う
4.1 COVID-19
患者の年代と公表曜日の関係を調べた。
library(ca)
library(readr)
df0 <- read.csv("./data/df0.csv")
df <- df0 # データの再読み込み
library(dplyr)
df$day <- recode(df$day,"月"="Mon","火"="Tue", "水"="Wed", "木"="Thu", "金"="Fri", "土"="Sat", "日"="Sun")
df$age <- recode(df$age,"'-"= NULL,
"100歳以上" = "a",
"10歳未満"="b",
"10代"="c",
"20代"="d",
"30代"="e",
"40代"="f",
"50代"="g",
"60代"="h",
"70代"="i",
"80代"="j",
"90代"="k",
"不明" = NULL)
df$sex <- recode(df$sex,"―"=NULL,"'-"=NULL,"女性"="F","男性"="M")
dat_ca <- table(df$age,df$day)
# a,b,c,kを除外する。
dat_ca <- dat_ca[-c(1:3,11),]
dat <- dat_ca
knitr::kable(dat)| Fri | Mon | Sat | Sun | Thu | Tue | Wed | |
|---|---|---|---|---|---|---|---|
| d | 676 | 410 | 624 | 535 | 752 | 475 | 482 |
| e | 453 | 237 | 401 | 344 | 519 | 319 | 321 |
| f | 262 | 168 | 288 | 208 | 297 | 239 | 222 |
| g | 216 | 133 | 254 | 146 | 204 | 216 | 168 |
| h | 125 | 70 | 133 | 107 | 138 | 118 | 114 |
| i | 106 | 58 | 116 | 104 | 127 | 93 | 99 |
| j | 91 | 41 | 65 | 72 | 75 | 52 | 57 |
# -------------
# 右側のオブジェクト名を変更する
# dat <- dat
# Package "ca"
# a scree-plot of the principal inertias
# and row and column contributions
summary(ca(dat))##
## Principal inertias (eigenvalues):
##
## dim value % cum% scree plot
## 1 0.004157 69.6 69.6 *****************
## 2 0.000914 15.3 84.9 ****
## 3 0.000580 9.7 94.6 **
## 4 0.000281 4.7 99.3 *
## 5 4e-05000 0.7 100.0
## 6 00000000 0.0 100.0
## -------- -----
## Total: 0.005973 100.0
##
##
## Rows:
## name mass qlt inr k=1 cor ctr k=2 cor ctr
## 1 | d | 343 905 128 | -40 706 130 | 21 199 166 |
## 2 | e | 225 721 122 | -48 721 126 | 0 0 0 |
## 3 | f | 146 864 76 | 52 864 94 | 0 0 0 |
## 4 | g | 116 972 401 | 139 939 541 | 26 33 86 |
## 5 | h | 70 969 73 | 52 440 46 | -57 529 251 |
## 6 | i | 61 907 76 | 5 4 0 | -82 904 447 |
## 7 | j | 39 408 125 | -81 347 62 | -34 61 50 |
##
## Columns:
## name mass qlt inr k=1 cor ctr k=2 cor ctr
## 1 | Fri | 167 561 112 | -43 471 76 | 19 90 66 |
## 2 | Mon | 97 743 79 | 2 1 0 | 60 743 384 |
## 3 | Sat | 163 972 138 | 69 952 189 | 10 20 18 |
## 4 | Sun | 131 889 151 | -66 641 139 | -41 248 244 |
## 5 | Thu | 183 735 200 | -69 726 208 | 8 9 12 |
## 6 | Tue | 131 985 266 | 109 979 374 | -9 6 11 |
## 7 | Wed | 127 927 55 | 22 190 15 | -44 737 264 |# using FactoMineR
# CA with function CA
library(FactoMineR)
# apply CA
ca = CA(dat, graph = T)
# matrix with eigenvalues
#ca$eig
# screeplot
draw_screeplot <- function(res.CA) {
library(qcc)
knitr::kable(res.CA$eig)
drange <- dim(res.CA$eig)[1]
dnames <- paste("Dim",1:drange,sep="")
pareto.chart(res.CA$eig[,1],names=dnames,
main="Screeplot",las=1,xlab="座標軸",
ylab="固有値",ylab2="積％",
col="lightblue")
} # Fugimoto's CA Tools
par(family="HiraKakuProN-W3")
draw_screeplot(ca)## Package 'qcc' version 2.7## Type 'citation("qcc")' for citing this R package in publications.
##
## Pareto chart analysis for res.CA$eig[, 1]
## Frequency Cum.Freq. Percentage Cum.Percent.
## dim 1 4.156830e-03 4.156830e-03 6.959946e+01 6.959946e+01
## dim 2 9.140420e-04 5.070872e-03 1.530417e+01 8.490363e+01
## dim 3 5.798275e-04 5.650700e-03 9.708283e+00 9.461192e+01
## dim 4 2.809882e-04 5.931688e-03 4.704697e+00 9.931661e+01
## dim 5 4.034623e-05 5.972034e-03 6.755330e-01 9.999215e+01
## dim 6 4.691003e-07 5.972503e-03 7.854333e-03 1.000000e+02# -----------
col <- ca$col$coord # 列スコア
colnames(col)[1:2] <- c("Dim_1","Dim_2") # 列名の修正
col <- as.data.frame(col[,c(1:2)]) # データフレーム化
col$c <- "1" # 色分け用
row <- ca$row$coord # 行スコア
colnames(row)[1:2] <- c("Dim_1","Dim_2") # 列名の修正
row <- as.data.frame(row[,c(1:2)])
row$c <- "2"
dat <- rbind(col,row) # 列スコアと行スコアの併合
code <- rownames(dat) # 表示用ラベル
library(patchwork)
library(ggrepel)
p1 <- ggplot(data = dat, aes(x = Dim_1, y = Dim_2,shape= c,colour=c)) +
geom_point()+
geom_hline(yintercept = 0, colour = "gray75",linetype="dashed") +
geom_vline(xintercept = 0, colour = "gray75",linetype="dashed") +
geom_text_repel(aes(label=code))+
scale_shape_manual(values=c(17,16))+
scale_colour_brewer(palette="Set1")+
# モノクロ用
#scale_colour_grey(start = 0,end=0.5)+
labs(x = "Dim 1", y = "Dim 2") +
ggtitle("CA plot (1)")+
theme_minimal()+
theme(legend.position = "none")
col <- ca$col$coord # 列スコア
colnames(col)[c(1,3)] <- c("Dim_1","Dim_3") # 列名の修正
col <- as.data.frame(col[,c(1,3)]) # データフレーム化
col$c <- "1" # 色分け用
row <- ca$row$coord # 行スコア
colnames(row)[c(1,3)] <- c("Dim_1","Dim_3") # 列名の修正
row <- as.data.frame(row[,c(1,3)])
row$c <- "2"
dat <- rbind(col,row) # 列スコアと行スコアの併合
code <- rownames(dat) # 表示用ラベル
p2 <- ggplot(data = dat, aes(x = Dim_1, y = Dim_3,shape= c,colour=c)) +
geom_point()+
geom_hline(yintercept = 0, colour = "gray75",linetype="dashed") +
geom_vline(xintercept = 0, colour = "gray75",linetype="dashed") +
geom_text_repel(aes(label=code))+
scale_shape_manual(values=c(17,16))+
scale_colour_brewer(palette="Set1")+
# モノクロ用
#scale_colour_grey(start = 0,end=0.5)+
labs(x = "Dim 1", y = "Dim 3") +
ggtitle("CA plot (2)")+
theme_minimal()+
theme(legend.position = "none")
p1
p2
4.1.1 symmetric, rowprincipal, colprincipal, symbiplot
dat <- dat_ca
par(mfrow=c(1,1))
plot(ca(dat),map="symmetric", main="symmetric")
plot(ca(dat),map="rowprincipal", main="rowprincipal:row=age")
plot(ca(dat),map="colprincipal",main="colprincipal:col=day")
plot(ca(dat),map="symbiplot",main="symbiplot")
cacoord(ca(dat),type= "symmetric",dim = c(1,2))## $rows
## Dim1 Dim2
## d -0.039704471 0.0210622628
## e -0.048306256 0.0004817867
## f 0.051779970 -0.0001111821
## g 0.139237738 0.0259917812
## h 0.052238253 -0.0573135409
## i 0.005251169 -0.0818852911
## j -0.081244558 -0.0339756344
##
## $columns
## Dim1 Dim2
## Fri -0.043323391 0.018957878
## Mon 0.001672622 0.060219704
## Sat 0.069334747 0.010130904
## Sun -0.066238742 -0.041217639
## Thu -0.068761559 0.007850891
## Tue 0.108830275 -0.008731660
## Wed 0.022129058 -0.043598263cacoord(ca(dat),type= "rowprincipal",dim = c(1,2))## $rows
## Dim1 Dim2
## d -0.039704471 0.0210622628
## e -0.048306256 0.0004817867
## f 0.051779970 -0.0001111821
## g 0.139237738 0.0259917812
## h 0.052238253 -0.0573135409
## i 0.005251169 -0.0818852911
## j -0.081244558 -0.0339756344
##
## $columns
## Dim1 Dim2
## Fri -0.67195672 0.6270565
## Mon 0.02594279 1.9918450
## Sat 1.07539942 0.3350928
## Sun -1.02737960 -1.3633270
## Thu -1.06650912 0.2596784
## Tue 1.68798502 -0.2888110
## Wed 0.34322728 -1.4420692cacoord(ca(dat),type= "colprincipal",dim = c(1,2))## $rows
## Dim1 Dim2
## d -0.61582636 0.696661716
## e -0.74924221 0.015935720
## f 0.80312039 -0.003677493
## g 2.15961245 0.859711946
## h 0.81022848 -1.895719857
## i 0.08144696 -2.708462429
## j -1.26012359 -1.123788264
##
## $columns
## Dim1 Dim2
## Fri -0.043323391 0.018957878
## Mon 0.001672622 0.060219704
## Sat 0.069334747 0.010130904
## Sun -0.066238742 -0.041217639
## Thu -0.068761559 0.007850891
## Tue 0.108830275 -0.008731660
## Wed 0.022129058 -0.043598263cacoord(ca(dat),type= "symbiplot",dim = c(1,2))## $rows
## Dim1 Dim2
## d -0.15636835 0.1211332825
## e -0.19024481 0.0027708514
## f 0.20392535 -0.0006394306
## g 0.54836079 0.1494839284
## h 0.20573021 -0.3296216280
## i 0.02068071 -0.4709386738
## j -0.31996591 -0.1954006632
##
## $columns
## Dim1 Dim2
## Fri -0.170620760 0.10903055
## Mon 0.006587296 0.34633555
## Sat 0.273061434 0.05826485
## Sun -0.260868421 -0.23705088
## Thu -0.270804043 0.04515204
## Tue 0.428606899 -0.05021752
## Wed 0.087150999 -0.250742324.1 2 矢印を入れる
par(pty="s")
plot(ca(dat_ca),map="symbiplot",main="symbiplot",arrows = c("TRUE","TRUE"),lwd=0.3)
par(pty="s")
plot(ca(dat_ca),map="symmetric",main="symmetric",arrows = c("TRUE","TRUE"),lwd=0.3)
4.2 Housing
- 『統計でみる日本2015』27ページの表
`
library(readr)
housing <- read_csv("./data/housing.csv")## New names:
## Rows: 4 Columns: 4
## ── Column specification
## ────────────────────────────────────────────────────────
## Delimiter: "," chr (1): ...1 dbl (3): A, B, C
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`housing <- housing[-1,-1]
housing <- data.frame(housing)
colnames(housing) <- c("単独世帯","核家族世帯","３世代世帯")
rownames(housing) <- c("一戸建","長屋建","共同住宅")
knitr::kable(housing)| 単独世帯 | 核家族世帯 | ３世代世帯 | |
|---|---|---|---|
| 一戸建 | 4674931 | 18822732 | 3413996 |
| 長屋建 | 517818 | 710774 | 27959 |
| 共同住宅 | 10827783 | 9564319 | 206872 |
housing## 単独世帯 核家族世帯 ３世代世帯
## 一戸建 4674931 18822732 3413996
## 長屋建 517818 710774 27959
## 共同住宅 10827783 9564319 206872library(ca)
ca(housing)##
## Principal inertias (eigenvalues):
## 1 2
## Value 0.157382 0.000209
## Percentage 99.87% 0.13%
##
##
## Rows:
## 一戸建 長屋建 共同住宅
## Mass 0.551840 0.025766 0.422394
## ChiDist 0.356133 0.244573 0.451380
## Inertia 0.069990 0.001541 0.086060
## Dim. 1 -0.897701 0.574651 1.137753
## Dim. 2 0.079080 -6.122101 0.270138
##
##
## Columns:
## 単独世帯 核家族世帯 ３世代世帯
## Mass 0.328511 0.596668 0.074821
## ChiDist 0.525647 0.192930 0.772183
## Inertia 0.090769 0.022209 0.044613
## Dim. 1 1.324855 -0.485719 -1.943509
## Dim. 2 0.537403 -0.663363 2.930519ca(t(housing))##
## Principal inertias (eigenvalues):
## 1 2
## Value 0.157382 0.000209
## Percentage 99.87% 0.13%
##
##
## Rows:
## 単独世帯 核家族世帯 ３世代世帯
## Mass 0.328511 0.596668 0.074821
## ChiDist 0.525647 0.192930 0.772183
## Inertia 0.090769 0.022209 0.044613
## Dim. 1 -1.324855 0.485719 1.943509
## Dim. 2 0.537403 -0.663363 2.930519
##
##
## Columns:
## 一戸建 長屋建 共同住宅
## Mass 0.551840 0.025766 0.422394
## ChiDist 0.356133 0.244573 0.451380
## Inertia 0.069990 0.001541 0.086060
## Dim. 1 0.897701 -0.574651 -1.137753
## Dim. 2 0.079080 -6.122101 0.270138cacoord: Extracting coordinates from ca objects.
cacoord(ca(housing),type= "symmetric",dim = c(1,2))## $rows
## Dim1 Dim2
## 一戸建 -0.3561311 0.001144043
## 長屋建 0.2279725 -0.088568087
## 共同住宅 0.4513632 0.003908069
##
## $columns
## Dim1 Dim2
## 単独世帯 0.5255892 0.007774572
## 核家族世帯 -0.1926916 -0.009596839
## ３世代世帯 -0.7710183 0.042395643par(family= "HiraKakuProN-W3")
par(mfrow=c(1,1),pty="s")
plot(ca(housing),main="列が家族類型") # 三角赤が家族類型、丸青が住居形態
plot(ca(t(housing)),main="行が家族類型") # 逆の設定
par(mfrow=c(1,2),pty="s",plt=c(0,1,0,0.65),family= "HiraKakuProN-W3")
mosaicplot(housing,color = TRUE,main = "住居形態 --> 家族形態")
mosaicplot(t(housing),color = TRUE,main = "家族形態 --> 住居形態")
par(mfrow=c(1,2),plt=c(0,1,0.1,0.9),pty="s",family= "HiraKakuProN-W3", cex=0.8) # 文字のサイズを指定(cex)
mosaicplot(t(housing),color = TRUE,main = "家族形態 --> 住居形態")
plot(ca(t(housing)), main = "対応分析")
5. 変数のカテゴリー化
# 分析用データの読み込み
dat <- read.csv("./data/dat_ssdse.csv")
head(dat)## X Code Prefecture Municipality A1101 A6108 B1103 A5101 A5102 D2201 F1102
## 1 1 1100 北海道 札幌市 1952356 100.4 43898 119314 110535 0.72 844313
## 2 2 1202 北海道 函館市 265979 102.8 12333 8748 9549 0.46 117125
## 3 3 1203 北海道 小樽市 121924 101.9 8009 3185 3751 0.42 51317
## 4 4 1204 北海道 旭川市 339605 100.6 35194 10245 11075 0.49 152385
## 5 5 1205 北海道 室蘭市 88564 109.4 4327 3061 3651 0.61 37286
## 6 6 1206 北海道 釧路市 174742 100.6 30910 5790 6813 0.44 74840
## F2701 A710201 A710202 A1303 A4101 A4200 A710101 A810101 A810102 A810103
## 1 372730 1899980 52376 483534 13821 18668 920415 531945 493644 38301
## 2 102578 255149 10830 85931 1410 3633 123651 74373 66690 7683
## 3 41516 116724 5200 45240 544 1885 55299 35155 31709 3446
## 4 139119 326243 13362 106444 2201 4186 155218 96178 88393 7785
## 5 31920 86029 2535 30118 515 1299 43536 24810 23014 1796
## 6 65204 169155 5587 52867 951 2194 81846 49477 45098 4379
## A810104 A810105 D2205 D2206 Dnsty In Out OutOverIn
## 1 12311 375242 239213005 326840745 44.474828 0.06111283 0.05661621 0.926421
## 2 999 48247 27485413 58883567 21.566448 0.03288981 0.03590133 1.091564
## 3 228 19911 11802308 27497809 15.223374 0.02612283 0.03076507 1.177708
## 4 1078 57488 34725118 67224590 9.649514 0.03016740 0.03261142 1.081015
## 5 239 18486 11239460 18688855 20.467761 0.03456258 0.04122443 1.192747
## 6 671 31697 18217586 40713325 5.653251 0.03313456 0.03898891 1.176684
## sc r_instHm r_nuclear r_nuclear_relativeH r_OnePerson aged
## 1 44.14595 2.682707 53.63276 0.9279982 0.4076878 24.76669
## 2 87.57994 4.071750 53.93406 0.8966964 0.3901869 32.30744
## 3 80.90107 4.264952 57.34100 0.9019770 0.3600608 37.10508
## 4 91.29442 3.934571 56.94765 0.9190563 0.3703694 31.34347
## 5 85.60854 2.862337 52.86200 0.9276098 0.4246141 34.00705
## 6 87.12453 3.197285 55.10104 0.9114942 0.3872761 30.25432
## naturalInc
## 1 -4847
## 2 -2223
## 3 -1341
## 4 -1985
## 5 -784
## 6 -12435.1 財政力指数
財政力指数と人口との関係(散布図)
library(ggrepel)
ggplot(dat,aes(x=A1101,y=D2201))+
geom_point()+
labs(x="人口(100万人以下)",y="財政力指数",title = "全国")+
scale_x_continuous(limits = c(0,1000000))+
theme_classic(base_family="HiraKakuProN-W3")## Warning: Removed 34 rows containing missing values (`geom_point()`).
人口を４レベルに分ける
pop <- as.numeric(dat$A1101)
# 21478.16 48434.88
#pop <- cut(pop,breaks = c(0,21478,48435,Inf),labels=c("small","medium","large"),include.lowest = TRUE)
#pop
#table(pop)
#table(cut(pop,quantile(pop),include.lowest = T))
#table(cut(pop,quantile(pop),include.lowest = F))
pop <- cut(pop,quantile(pop),labels=c("p1","p2","p3","p4"),include.lowest = TRUE)
#pop
table(pop)## pop
## p1 p2 p3 p4
## 436 435 435 435財政力指数を４レベルに分ける
fscl <- as.numeric(dat$D2201)
#table(cut(fscl,quantile(fscl),include.lowest = T))
#table(cut(fscl,quantile(fscl),include.lowest = F))
fscl <- cut(fscl,quantile(fscl,na.rm = TRUE),labels=c("f1","f2","f3","f4"),include.lowest = TRUE,na.rm =TRUE)
#fscl
table(fscl)## fscl
## f1 f2 f3 f4
## 441 423 427 427table(fscl,pop)## pop
## fscl p1 p2 p3 p4
## f1 328 94 19 0
## f2 73 188 139 23
## f3 14 109 173 131
## f4 21 44 103 259#tbl <- table(fscl,pop)
#write.csv(tbl,file = "tbl.csv")
#
#test <- read.csv("tbl.csv",row.names = 1)
#test
library(vcd)## 要求されたパッケージ grid をロード中ですmosaic(table(pop,fscl),shade=T)
par(family="HiraKakuProN-W3")
mosaicplot(table(pop,fscl),shade=T,main="人口と財政力指数",sub = "p1 < p2 < p3 < p4")
pop <- as.numeric(dat$A1101)
# 21478.16 48434.88
#pop <- cut(pop,breaks = c(0,21478,48435,Inf),labels=c("small","medium","large"),include.lowest = TRUE)
#pop
#table(pop)
#table(cut(pop,quantile(pop),include.lowest = T))
#table(cut(pop,quantile(pop),include.lowest = F))
pop <- cut(pop,quantile(pop),labels=c("p1","p2","p3","p4"),include.lowest = TRUE)
#pop
table(pop)## pop
## p1 p2 p3 p4
## 436 435 435 435fscl <- as.numeric(dat$D2201)
#table(cut(fscl,quantile(fscl),include.lowest = T))
#table(cut(fscl,quantile(fscl),include.lowest = F))
fscl <- cut(fscl,quantile(fscl,na.rm = TRUE),labels=c("f1","f2","f3","f4"),include.lowest = TRUE,na.rm =TRUE)
#fscl
table(fscl)## fscl
## f1 f2 f3 f4
## 441 423 427 427table(fscl,pop)## pop
## fscl p1 p2 p3 p4
## f1 328 94 19 0
## f2 73 188 139 23
## f3 14 109 173 131
## f4 21 44 103 259#tbl <- table(fscl,pop)
#write.csv(tbl,file = "tbl.csv")
#
#test <- read.csv("tbl.csv",row.names = 1)
#test
library(vcd)
mosaic(table(pop,fscl),shade=T)
par(family="HiraKakuProN-W3")
mosaicplot(table(pop,fscl),shade=T,main="人口と財政力指数",sub = "p1 < p2 < p3 < p4")
5.2 転出超過率
OutOverIn <- as.numeric(dat$OutOverIn)
OutOverIn <- cut(OutOverIn,quantile(OutOverIn,na.rm = TRUE),labels=c("o1","o2","o3","o4"),include.lowest = TRUE,na.rm =TRUE)
#OutOverIn
table(OutOverIn)## OutOverIn
## o1 o2 o3 o4
## 436 435 435 435table(OutOverIn,pop)## pop
## OutOverIn p1 p2 p3 p4
## o1 89 64 94 189
## o2 69 88 126 152
## o3 90 126 135 84
## o4 188 157 80 10summary(table(OutOverIn,pop))## Number of cases in table: 1741
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 315.82, df = 9, p-value = 1.146e-62library(vcd)
mosaic(table(pop,OutOverIn),shade=T)
par(family="HiraKakuProN-W3")
mosaicplot(table(pop,OutOverIn),shade=T,main="人口と転出超過率",sub = "p1 < p2 < p3 < p4")
5.3 核家族率
5.3.1 モザイク図
r_nuclear_relativeH <- as.numeric(dat$r_nuclear_relativeH)
r_nuclear_relativeH <- cut(r_nuclear_relativeH,quantile(r_nuclear_relativeH,na.rm = TRUE),labels=c("o1","o2","o3","o4"),include.lowest = TRUE,na.rm =TRUE)
#r_nuclear_relativeH
table(r_nuclear_relativeH)## r_nuclear_relativeH
## o1 o2 o3 o4
## 434 434 434 434table(r_nuclear_relativeH,pop)## pop
## r_nuclear_relativeH p1 p2 p3 p4
## o1 139 162 99 34
## o2 103 118 140 73
## o3 120 97 110 107
## o4 69 58 86 221summary(table(r_nuclear_relativeH,pop))## Number of cases in table: 1736
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 270.07, df = 9, p-value = 5.732e-53library(vcd)
mosaic(table(pop,r_nuclear_relativeH),shade=T)
par(family="HiraKakuProN-W3")
mosaicplot(table(pop,r_nuclear_relativeH),shade=T,main="人口と核家族率",sub = "p1 < p2 < p3 < p4")
5.3.2 三重クロス集計
table(OutOverIn,r_nuclear_relativeH,pop)## , , pop = p1
##
## r_nuclear_relativeH
## OutOverIn o1 o2 o3 o4
## o1 22 16 28 23
## o2 22 13 21 13
## o3 24 25 24 17
## o4 71 49 47 16
##
## , , pop = p2
##
## r_nuclear_relativeH
## OutOverIn o1 o2 o3 o4
## o1 12 13 28 11
## o2 23 36 20 9
## o3 49 33 26 18
## o4 78 36 23 20
##
## , , pop = p3
##
## r_nuclear_relativeH
## OutOverIn o1 o2 o3 o4
## o1 1 27 36 30
## o2 16 41 35 34
## o3 40 47 30 18
## o4 42 25 9 4
##
## , , pop = p4
##
## r_nuclear_relativeH
## OutOverIn o1 o2 o3 o4
## o1 3 15 32 139
## o2 7 36 50 59
## o3 23 18 22 21
## o4 1 4 3 2summary(table(OutOverIn,r_nuclear_relativeH,pop))## Number of cases in table: 1736
## Number of factors: 3
## Test for independence of all factors:
## Chisq = 1068, df = 54, p-value = 2.57e-188library(vcd)
mosaic(table(pop,r_nuclear_relativeH,OutOverIn),shade=T)
par(family="HiraKakuProN-W3")
mosaic(table(pop,r_nuclear_relativeH,OutOverIn),shade=T,direction = c("v","v","h"))
mosaic(table(pop,r_nuclear_relativeH,OutOverIn),shade=T,direction = c("v","h","v"))
mosaic(table(pop,r_nuclear_relativeH,OutOverIn),shade=T,direction = c("h","v","v"))
mosaic(table(pop,r_nuclear_relativeH,OutOverIn),shade=T,direction = c("v","h","h"))
mosaic(table(pop,r_nuclear_relativeH,OutOverIn),shade=T,direction = c("h","v","h"))
mosaic(table(pop,r_nuclear_relativeH,OutOverIn),shade=T,direction = c("h","h","v"))
par(family="HiraKakuProN-W3")
mosaicplot(table(pop,r_nuclear_relativeH,OutOverIn),shade=T,main="転出超過率",sub = "")
5.4 自地域就業者率
5.4.1 モザイク図
sc <- as.numeric(dat$sc)
sc <- cut(sc,quantile(sc,na.rm = TRUE),labels=c("s1","s2","s3","s4"),include.lowest = TRUE,na.rm =TRUE)
#sc
table(sc)## sc
## s1 s2 s3 s4
## 434 434 434 434table(sc,pop)## pop
## sc p1 p2 p3 p4
## s1 43 112 120 159
## s2 84 132 125 93
## s3 119 111 114 90
## s4 185 80 76 93summary(table(sc,pop))## Number of cases in table: 1736
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 158.14, df = 9, p-value = 1.808e-29library(vcd)
mosaic(table(pop,sc),shade=T)
par(family="HiraKakuProN-W3")
mosaicplot(table(pop,sc),shade=T,main="自地域就業者率",sub = "p1 < p2 < p3 < p4")
library(ggrepel)
ggplot(dat,aes(x=A1101,y=sc))+
geom_point()+
labs(x="人口(100万人以下)",y="自地域就業者率",title = "全国")+
scale_x_continuous(limits = c(0,1000000))+
theme_classic(base_family="HiraKakuProN-W3")## Warning: Removed 16 rows containing missing values (`geom_point()`).
5.4.2 対応分析
tbl <- table(sc,pop)
library(ca)
ca(tbl)##
## Principal inertias (eigenvalues):
## 1 2 3
## Value 0.077691 0.013251 0.000152
## Percentage 85.29% 14.55% 0.17%
##
##
## Rows:
## s1 s2 s3 s4
## Mass 0.250000 0.250000 0.250000 0.250000
## ChiDist 0.382178 0.185380 0.104146 0.416059
## Inertia 0.036515 0.008591 0.002712 0.043276
## Dim. 1 -1.284668 -0.410692 0.237828 1.457532
## Dim. 2 1.160296 -1.261061 -0.677065 0.777830
##
##
## Columns:
## p1 p2 p3 p4
## Mass 0.248272 0.250576 0.250576 0.250576
## ChiDist 0.483373 0.170967 0.177518 0.267013
## Inertia 0.058009 0.007324 0.007896 0.017865
## Dim. 1 1.733128 -0.454380 -0.557632 -0.705179
## Dim. 2 0.146962 -0.987213 -0.728256 1.569859plot(ca(tbl))
library(FactoMineR)
CA(tbl)
## **Results of the Correspondence Analysis (CA)**
## The row variable has 4 categories; the column variable has 4 categories
## The chi square of independence between the two variables is equal to 158.1398 (p-value = 1.807913e-29 ).
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$col" "results for the columns"
## 3 "$col$coord" "coord. for the columns"
## 4 "$col$cos2" "cos2 for the columns"
## 5 "$col$contrib" "contributions of the columns"
## 6 "$row" "results for the rows"
## 7 "$row$coord" "coord. for the rows"
## 8 "$row$cos2" "cos2 for the rows"
## 9 "$row$contrib" "contributions of the rows"
## 10 "$call" "summary called parameters"
## 11 "$call$marge.col" "weights of the columns"
## 12 "$call$marge.row" "weights of the rows"5.5 カテゴリー数を増やす
各変数を4分割ではなく5分割して、モザイク図を描く。また、対応分析をおこなう。ある程度の数がある場合には、その方がよい。
pop <- as.numeric(dat$A1101)
qa <- quantile(pop,c(0,0.20,0.40,0.60,0.80,1),na.rm = TRUE)
pop5 <- cut(pop,qa,labels = c("p1","p2","p3","p4","p5"),include.lowest = TRUE)
table(pop5)## pop5
## p1 p2 p3 p4 p5
## 349 348 348 348 348sc <- as.numeric(dat$sc)
qa <- quantile(sc,c(0,0.20,0.40,0.60,0.80,1),na.rm = TRUE)
sc5 <- cut(sc,qa,labels = c("s1","s2","s3","s4","s5"),include.lowest = TRUE)
table(sc5)## sc5
## s1 s2 s3 s4 s5
## 348 347 347 347 347table(sc5,pop5)## pop5
## sc5 p1 p2 p3 p4 p5
## s1 16 61 78 85 108
## s2 38 95 82 73 59
## s3 65 72 77 80 53
## s4 91 62 66 61 67
## s5 134 58 45 49 61summary(table(sc5,pop5))## Number of cases in table: 1736
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 189.66, df = 16, p-value = 9.657e-32library(vcd)
mosaic(table(pop5,sc5),shade=T)
par(family="HiraKakuProN-W3")
mosaicplot(table(pop5,sc5),shade=T,main="自地域就業者率",sub = "5分類")
tbl <- table(sc5,pop5)
library(ca)
ca(tbl)##
## Principal inertias (eigenvalues):
## 1 2 3 4
## Value 0.089996 0.016827 0.002172 0.000256
## Percentage 82.38% 15.4% 1.99% 0.23%
##
##
## Rows:
## s1 s2 s3 s4 s5
## Mass 0.200461 0.199885 0.199885 0.199885 0.199885
## ChiDist 0.438895 0.279270 0.137875 0.164102 0.479018
## Inertia 0.038614 0.015589 0.003800 0.005383 0.045865
## Dim. 1 -1.319432 -0.690975 -0.092164 0.527535 1.578838
## Dim. 2 1.461076 -1.366727 -0.830601 0.231477 0.500563
##
##
## Columns:
## p1 p2 p3 p4 p5
## Mass 0.198157 0.200461 0.200461 0.200461 0.200461
## ChiDist 0.599733 0.194999 0.192141 0.187210 0.281898
## Inertia 0.071273 0.007622 0.007401 0.007026 0.015930
## Dim. 1 1.998213 -0.272841 -0.582457 -0.578779 -0.541169
## Dim. 2 0.125328 -1.265627 -0.537444 -0.074874 1.754057plot(ca(tbl))
library(FactoMineR)
CA(tbl)
## **Results of the Correspondence Analysis (CA)**
## The row variable has 5 categories; the column variable has 5 categories
## The chi square of independence between the two variables is equal to 189.6607 (p-value = 9.656683e-32 ).
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$col" "results for the columns"
## 3 "$col$coord" "coord. for the columns"
## 4 "$col$cos2" "cos2 for the columns"
## 5 "$col$contrib" "contributions of the columns"
## 6 "$row" "results for the rows"
## 7 "$row$coord" "coord. for the rows"
## 8 "$row$cos2" "cos2 for the rows"
## 9 "$row$contrib" "contributions of the rows"
## 10 "$call" "summary called parameters"
## 11 "$call$marge.col" "weights of the columns"
## 12 "$call$marge.row" "weights of the rows"5.6 応用(流出超過率)
# 分析用データの読み込み
dat <- read.csv("./data/dat_ssdse.csv")
head(dat)## X Code Prefecture Municipality A1101 A6108 B1103 A5101 A5102 D2201 F1102
## 1 1 1100 北海道 札幌市 1952356 100.4 43898 119314 110535 0.72 844313
## 2 2 1202 北海道 函館市 265979 102.8 12333 8748 9549 0.46 117125
## 3 3 1203 北海道 小樽市 121924 101.9 8009 3185 3751 0.42 51317
## 4 4 1204 北海道 旭川市 339605 100.6 35194 10245 11075 0.49 152385
## 5 5 1205 北海道 室蘭市 88564 109.4 4327 3061 3651 0.61 37286
## 6 6 1206 北海道 釧路市 174742 100.6 30910 5790 6813 0.44 74840
## F2701 A710201 A710202 A1303 A4101 A4200 A710101 A810101 A810102 A810103
## 1 372730 1899980 52376 483534 13821 18668 920415 531945 493644 38301
## 2 102578 255149 10830 85931 1410 3633 123651 74373 66690 7683
## 3 41516 116724 5200 45240 544 1885 55299 35155 31709 3446
## 4 139119 326243 13362 106444 2201 4186 155218 96178 88393 7785
## 5 31920 86029 2535 30118 515 1299 43536 24810 23014 1796
## 6 65204 169155 5587 52867 951 2194 81846 49477 45098 4379
## A810104 A810105 D2205 D2206 Dnsty In Out OutOverIn
## 1 12311 375242 239213005 326840745 44.474828 0.06111283 0.05661621 0.926421
## 2 999 48247 27485413 58883567 21.566448 0.03288981 0.03590133 1.091564
## 3 228 19911 11802308 27497809 15.223374 0.02612283 0.03076507 1.177708
## 4 1078 57488 34725118 67224590 9.649514 0.03016740 0.03261142 1.081015
## 5 239 18486 11239460 18688855 20.467761 0.03456258 0.04122443 1.192747
## 6 671 31697 18217586 40713325 5.653251 0.03313456 0.03898891 1.176684
## sc r_instHm r_nuclear r_nuclear_relativeH r_OnePerson aged
## 1 44.14595 2.682707 53.63276 0.9279982 0.4076878 24.76669
## 2 87.57994 4.071750 53.93406 0.8966964 0.3901869 32.30744
## 3 80.90107 4.264952 57.34100 0.9019770 0.3600608 37.10508
## 4 91.29442 3.934571 56.94765 0.9190563 0.3703694 31.34347
## 5 85.60854 2.862337 52.86200 0.9276098 0.4246141 34.00705
## 6 87.12453 3.197285 55.10104 0.9114942 0.3872761 30.25432
## naturalInc
## 1 -4847
## 2 -2223
## 3 -1341
## 4 -1985
## 5 -784
## 6 -1243pop <- as.numeric(dat$A1101)
qa <- quantile(pop,c(0,0.20,0.40,0.60,0.80,1),na.rm = TRUE)
qa## 0% 20% 40% 60% 80% 100%
## 0 6505 15938 35206 80826 3724844# pop5 <- cut(pop,qa,labels = c("p1","p2","p3","p4","p5"),include.lowest = TRUE)
pop5 <- cut(pop,qa,labels = c("p1","p2","p3","p4","p5"),right=TRUE,include.lowest = TRUE)
table(pop5)## pop5
## p1 p2 p3 p4 p5
## 349 348 348 348 348OutOverIn <- as.numeric(dat$OutOverIn)
OOI5 <- cut(OutOverIn,quantile(OutOverIn,c(0,0.20,0.40,0.60,0.80,1),na.rm = TRUE),labels=c("OOI1","OOI2","OOI3","OOI4","OOI5"),include.lowest = TRUE,na.rm =TRUE)
#OutOverIn
table(OOI5)## OOI5
## OOI1 OOI2 OOI3 OOI4 OOI5
## 349 348 348 348 348table(OOI5,pop5)## pop5
## OOI5 p1 p2 p3 p4 p5
## OOI1 60 38 51 76 124
## OOI2 42 48 58 76 124
## OOI3 55 70 76 77 70
## OOI4 53 78 93 94 30
## OOI5 139 114 70 25 0summary(table(OOI5,pop5))## Number of cases in table: 1741
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 370.3, df = 16, p-value = 6.029e-69library(vcd)
mosaic(table(pop5,OOI5),shade=T)
par(family="HiraKakuProN-W3")
mosaicplot(table(pop5,OOI5),shade=T,main="人口規模別流出超過率",sub = "人口5分類")
tbl <- table(OOI5,pop5)
library(ca)
ca(tbl)##
## Principal inertias (eigenvalues):
## 1 2 3 4
## Value 0.17455 0.03665 0.001483 3e-06
## Percentage 82.07% 17.23% 0.7% 0%
##
##
## Rows:
## OOI1 OOI2 OOI3 OOI4 OOI5
## Mass 0.200460 0.199885 0.199885 0.199885 0.199885
## ChiDist 0.427070 0.424985 0.113973 0.355720 0.748991
## Inertia 0.036562 0.036102 0.002596 0.025293 0.112133
## Dim. 1 -0.935062 -0.994169 -0.127221 0.321954 1.737184
## Dim. 2 -0.848435 -0.396548 0.498640 1.715050 -0.966268
##
##
## Columns:
## p1 p2 p3 p4 p5
## Mass 0.200460 0.199885 0.199885 0.199885 0.199885
## ChiDist 0.503109 0.381167 0.210799 0.335120 0.712856
## Inertia 0.050740 0.029041 0.008882 0.022448 0.101574
## Dim. 1 1.053954 0.900971 0.251225 -0.568976 -1.640203
## Dim. 2 -1.250499 0.108812 0.952957 1.213339 -1.021016plot(ca(tbl))
library(FactoMineR)
CA(tbl)
## **Results of the Correspondence Analysis (CA)**
## The row variable has 5 categories; the column variable has 5 categories
## The chi square of independence between the two variables is equal to 370.2855 (p-value = 6.029486e-69 ).
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$col" "results for the columns"
## 3 "$col$coord" "coord. for the columns"
## 4 "$col$cos2" "cos2 for the columns"
## 5 "$col$contrib" "contributions of the columns"
## 6 "$row" "results for the rows"
## 7 "$row$coord" "coord. for the rows"
## 8 "$row$cos2" "cos2 for the rows"
## 9 "$row$contrib" "contributions of the rows"
## 10 "$call" "summary called parameters"
## 11 "$call$marge.col" "weights of the columns"
## 12 "$call$marge.row" "weights of the rows"6. べき乗分布への当てはめ
(別のファイルからのコピー。ファイル名は、日本の人口である。)
全国市町村の人口分布にべき乗分布、対数正規分布を当てはめる。
まずデータの読み込む。（読み込むファイルがどこにあるかを確認すること。また、working directoryに注意。「..」は、1つ上位のディレクトリを意味する。 ）
6.1 準備
# 分析用データの読み込み
dat <- read.csv("./data/dat_ssdse.csv")arrannge関数を使って、人口の少ない順にデータフレーム全体を並び替える。結果をdat_oとする。
library(dplyr)
# 人口の少ない順に並び替え
dat_o <- arrange(dat,A1101)ここからは、人口の列だけを取り出して処理する。 rev関数を使って人口(A1101)を多い順にならべる。popというオブジェクトとして保存する。
# 人口(A1101)のみを取り出す
rv_dt <- rev(dat_o$A1101) # 多い順に並べた人口(A1101)
pop <- rv_dt人口の順位をrankとする。rankに1から1741までの数が入る。
# ランク
rank <- 1:1741 # 人口の順位人口順位と人口のそれぞれを対数変換したデータに回帰分析をおこなって、両対数グラフにプロットしたときに直線となる関係を想定してその回帰直線の係数を求める。その係数からべき乗分布の計算式を導き出す。べき乗分布は、両対数グラフに描くと直線となるという特徴がある。
べき乗分布の計算は以下のようにしておこなった。最初のものは、上位100位までの人口に焦点を合わせ、2番目は、上位20位までに焦点を合わせた。
以下のものが計算されたべき乗分布の式である。
- y = (10^6.64)*x^(-0.6013)
- y = (10^6.6162)*x^(-0.5711)
x <- log10(rank[1:100])
y <- log10(pop[1:100])
df1 <- data.frame(x,y)
lm(y~x,data=df1)
Call:
lm(formula = y ~ x, data = df1)
Coefficients:
(Intercept) x
6.6430 -0.6013
x <- log10(rank[1:20])
y <- log10(pop[1:20])
df1 <- data.frame(x,y)
lm(y~x,data=df1)
Call:
lm(formula = y ~ x, data = df1)
Coefficients:
(Intercept) x
6.6162 -0.5711 # --------------------------
# べき乗分布の計算式(1)
e <- function(x) {
return(y <- (10^6.64)*x^(-0.6013))
}
y <- 1:1741
for (i in 1:1741) {
y[i] <- e(i)
}
# 計算式から出た数値
pop_e <- y# --------------------------
# べき乗分布の計算式(2)
e <- function(x) {
return(y <- (10^6.6162)*x^(-0.5711))
}
y <- 1:1741
for (i in 1:1741) {
y[i] <- e(i)
}
# 計算式から出た数値
pop_e_20 <- y
# --------------------------以下の5変数を集めたデータフレーム(df)をつくる。
x <- rank 人口のランク
y1 <- rv_dt 人口
y2 <- pop_e べき乗分布の推定値（100位までに最適化）
y3 <- pop_e_20 同上（20位までに最適化）
次のy4は、対数正規分布を想定した数値
y4 <- qlnorm(ppoints(1741),6,0.864,lower.tail = F)*480x <- rank
y1 <- rv_dt
y2 <- pop_e
y3 <- pop_e_20
y4 <- qlnorm(ppoints(1741),6,0.864,lower.tail = F)*480
df <- data.frame(x,y1,y2,y3,y4)
head(df)## x y1 y2 y3 y4
## 1 1 3724844 4365158 4132378 3793921
## 2 2 2691185 2877337 2781517 2904430
## 3 3 2295638 2254795 2206562 2544586
## 4 4 1952356 1896624 1872249 2324725
## 5 5 1538681 1658477 1648232 2168992
## 6 6 1537272 1486270 1485244 2049674tail(df)## x y1 y2 y3 y4
## 1736 1736 41 49210.61 58357.03 18294.96
## 1737 1737 18 49193.57 58337.84 17288.54
## 1738 1738 0 49176.55 58318.67 16130.38
## 1739 1739 0 49159.55 58299.51 14736.66
## 1740 1740 0 49142.56 58280.37 12910.87
## 1741 1741 0 49125.58 58261.25 9883.89long formatに変換しdf_longと名付ける。
df_long <- gather(df,"bucket","y",-x)
head(df_long)## x bucket y
## 1 1 y1 3724844
## 2 2 y1 2691185
## 3 3 y1 2295638
## 4 4 y1 1952356
## 5 5 y1 1538681
## 6 6 y1 1537272tail(df_long)## x bucket y
## 6959 1736 y4 18294.96
## 6960 1737 y4 17288.54
## 6961 1738 y4 16130.38
## 6962 1739 y4 14736.66
## 6963 1740 y4 12910.87
## 6964 1741 y4 9883.896.2 グラフ
べき乗分布、対数正規分布への当てはめ。 ggplot2でグラフを描く。
7. 対数正規分布への当てはめ
対数正規分布は、その分布を示す変数の値の対数を取ったものが正規分布を示すものである。対数を取ることができるのは正の数だけであるので、対数正規分布も、正の数だけを取ることができる。
7.1 準備
全国の市町村の人口について、対数正規分布を仮定しqlnorm関数を使い100個の推測値を発生させる。 右端の数値1000は、任意のもの。結果の数値を1000倍にしている。
y2 <- qlnorm(ppoints(100),0.1,3.155,lower.tail = F)*1000
y3 <- qlnorm(ppoints(100),6,0.864,lower.tail = F)*1000## x y1 y2 y3
## 1 1 3724844 3739628.1 3735050
## 2 2 2691185 1039645.6 2630576
## 3 3 2295638 535760.0 2193840
## 4 4 1952356 335821.0 1930416
## 5 5 1538681 232520.8 1745548
## 6 6 1537272 171108.8 1604936y1は実際の人口。それに近づくようにqlnorm関数のmeanlogとsdlogの値を調整してみた。1位の人口3724844に合わせた。
7.2 グラフを描く。
（include=TRUE,echo=FALSEと指定しているので結果だけが表示される。）
次のグラフでは、以下の2行が挿入され、x及びyの両軸が対数軸となっている。
scale_x_log10()
scale_y_log10()ggplot(df_long,aes(x,y,col=key))+
geom_point()+
geom_line()+
labs(title = "",subtitle = "y: 人口\ny1: 実際値\ny2: qlnorm(ppoints(100),0.1,3.155,lower.tail = F)*1000\ny3: qlnorm(ppoints(100),6,0.864,lower.tail = F)*1000 ", x="x: ランク",y="y",colour="凡例")+
scale_x_log10()+
scale_y_log10()+
theme_minimal(base_family="HiraKakuProN-W3")
次のグラフは、x軸とy軸とを関数coord_flipで反転させている。
以上のとおり、対数正規分布としては、y3が最も実際のデータと近い。
7.3 式の検討
y3 <- qlnorm(ppoints(100),6,0.864,lower.tail = F)*1000 ここでは、qlnormという関数を使っている。この関数は、対数正規分布(log-normal distribution)において、 「確率ベクトル」に対してそれにあたる値(分布曲線のx軸上の値)を計算するものである。 得られる数値は、quantile(分位点あるいは、分位数)とよばれる。 正規分布でこれに対応する関数はqnormである。
正規分布において、qnorm関数がどういう数値を計算するかは以下の以下の通りである。統計学でよく出てくる数値を確認することができる。
qnorm(0.05) = -1.644854
qnorm(0.05,lower.tail = F) = 1.644854qnormと対になっている関数は、pnormである。（対数正規分布では、plnorm。）
pnorm関数は、累積確率を計算するものである。累積確率は、lower.tail=FAlSEと指定することで上位確率を指定することになる。
以下のように、正規分布で中心は累積確率が0.5（つまり50パーセント）の点であり、上位確率5パーセント点は、1.64である。
pnorm(0)= 0.5
pnorm(1.644854) = 0.95
pnorm(1.644854,lower.tail = F) = 0.04999996
pnorm(-1.644854) = 0.04999996
pnorm(-1.644854,lower.tail = F) = 0.95人口の分布を正規分布などに当てはめる場合、qnormとpnormのどちらを使うことになるのか。 ランクから人口を推定するとしたら、ランクは上位から数えた累積確率とみることができるので、正規分布であればqnormを使えば、上位から5パーセントのところが測定値の1.64にあたると考え横軸に順位、縦軸に測定値の推定値を描くことができる。
正規分布の場合と同じように、
対数正規分布であればqnormに引数としては、ゼロから１までの数を累積確率として代入すればよい。
確率は大きい順に代入するために、rev(ppoints(100))とする必要はない。 lower.tailをFとすれば、順番に上位確率が想定されるので、そうすべきではない。
対数正規分布を当てはめる場合、実際の分布に合わせてmeanlogやsdlogの値を調整することになる。 最大値を実際の値に簡便に合わせたりする場合には、定数を乗じることも可能である。（対数正規分布では実際にゼロの値をとることはないが、原点のある「比率尺度」と考えていいだろう。）
7.4 meanlogとsdlogの意味
meanlogとsdlogを変化させると対数正規分布がどのように変化するかを調べてみてみよう。
y1とy3は、meanlogが同じ。y1とy2は、sdlogが同じ。 ここでは、y1は、これまでの当てはめで、最もよい式をつかっている。
y1 <- qlnorm(x,6,0.864,lower.tail = F)*1000これを基準として、y2は、meanlogを１増やしたもの、y3は、sdlogを1.2に増やしたものである。
