必要なパッケージ

require(tidyverse)
##  要求されたパッケージ tidyverse をロード中です
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.8     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.0
## ✔ readr   2.1.2     ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
require(maps)
##  要求されたパッケージ maps をロード中です 
## 
##  次のパッケージを付け加えます: 'maps' 
## 
##  以下のオブジェクトは 'package:purrr' からマスクされています: 
## 
##      map
require(plotly)
##  要求されたパッケージ plotly をロード中です 
## 
##  次のパッケージを付け加えます: 'plotly' 
## 
##  以下のオブジェクトは 'package:ggplot2' からマスクされています: 
## 
##      last_plot 
## 
##  以下のオブジェクトは 'package:stats' からマスクされています: 
## 
##      filter 
## 
##  以下のオブジェクトは 'package:graphics' からマスクされています: 
## 
##      layout
require(geosphere)
##  要求されたパッケージ geosphere をロード中です

データの準備

データを読み込む。

capital <- read.csv("prefCapital.csv")

中身を確認する。

head(capital)
##   Prefecture  Capital Longitude Latitude
## 1   Hokkaido  Sapporo  141.3469 43.06417
## 2     Aomori   Aomori  140.7400 40.82444
## 3      Iwate  Morioka  141.1525 39.70361
## 4     Miyagi   Sendai  140.8719 38.26889
## 5      Akita    Akita  140.1025 39.71861
## 6   Yamagata Yamagata  140.3633 38.24056

47都道府県の県庁所在地の情報である。

world <- map_data("world")

データを図示する。

world %>%
  ggplot(aes(x = long, y = lat, group = group)) + 
  geom_polygon(fill = "lightgrey", colour = "black", size = 0.1)

日本だけ抜き出してくる。

japan <- world %>% filter(region == "Japan")

日本地図を描く。

japan %>%
  ggplot(aes(x = long, y = lat, group =  group)) +
  geom_polygon(fill = "lightgray", colour = "black", size = 0.1)

日本地図を描く。

japan %>% 
  ggplot(aes(x = long, y = lat, group = group)) +
  geom_polygon(fill = "lightgray", colour = "black", size = 0.1) +
  geom_point(data = capital, mapping = aes(x = Longitude, y = Latitude, group = Prefecture))

##　主成分分析

主成分分析を行うデータ（県庁所在地の緯度経度）を準備する。

geo <- capital %>% select(Capital, Longitude, Latitude) %>% column_to_rownames(var = "Capital")
geo <- as.matrix(geo)
head(geo)
##          Longitude Latitude
## Sapporo   141.3469 43.06417
## Aomori    140.7400 40.82444
## Morioka   141.1525 39.70361
## Sendai    140.8719 38.26889
## Akita     140.1025 39.71861
## Yamagata  140.3633 38.24056

主成分分析を行う。

pca <- prcomp(geo)
summary(pca)
## Importance of components:
##                           PC1     PC2
## Standard deviation     4.3639 1.22318
## Proportion of Variance 0.9272 0.07284
## Cumulative Proportion  0.9272 1.00000

第1主成分の寄与率が93%、第2主成分は7％である。

主成分得点を図示する。

df <- data.frame(Capital = rownames(geo), pca$x)
plot_ly(data = df, x = ~PC1, y = ~PC2, text = ~Capital,
        type = "scatter", mode = "markers")

バイプロットを描く。

biplot(pca)

第1主成分は、緯度経度の両方が小さく（大きく）なると得点が大きく（小さく）なる。

第2主成分は、経度が大きく（小さく）緯度が小さい（大きい）と得点が大きく（小さく）なる。

固有値と固有ベクトルを確認してみる。

pca$sdev^2
## [1] 19.043569  1.496168
pca$rotation
##                  PC1        PC2
## Longitude -0.8333102  0.5528057
## Latitude  -0.5528057 -0.8333102

この計算は、以下のようにして行うこともできる。

まずは、分散共分散行列を求める。

V <- var(geo) # 分散共分散行列を計算する
V
##           Longitude Latitude
## Longitude 13.681186 8.083362
## Latitude   8.083362 6.858551

次に、分散共分散行列の固有値分解を行う。

eig <- eigen(V) # 固有値分解をする
eig # 結果を表示する
## eigen() decomposition
## $values
## [1] 19.043569  1.496168
## 
## $vectors
##            [,1]       [,2]
## [1,] -0.8333102  0.5528057
## [2,] -0.5528057 -0.8333102

主成分スコアを計算する。スコアを計算する前に、各変数（緯度・経度）の平均を0に基準化する。

X <- scale(geo, center = T, scale = F)
pcs <- X %*% eig$vectors # 主成分得点の計算
head(pcs)
##               [,1]         [,2]
## Sapporo  -8.668767 -3.467641518
## Aomori   -6.924862 -1.936771640
## Morioka  -6.649001 -0.774740245
## Sendai   -5.622086  0.265731346
## Akita    -5.782318 -1.367685904
## Yamagata -5.182596  0.008176506

主成分得点を図示する。

plot_ly(x = pcs[,1], y = pcs[,2], text = rownames(pcs), 
        type = "scatter", mode = "markers")

prcompを使って計算した場合と同じ結果となる（場合によっては、正負が逆転する）。

古典的多次元尺度法

多次元尺度法を適用する前に都市間の距離行列を求める。

緯度経度で表された2点間の距離は、geosphereパッケージのdistGeo関数で求められる。

distGeo(p1 = geo["Tokyo", ], p2 = geo["Osaka", ])
## [1] 395874.3

結果はメートルで表され、東京-大阪間はおよそ396 kmあることがわかる。

次に、全ての都市間で距離行列を計算する。

dm <- distm(geo, fun = distGeo)
rownames(dm) <- colnames(dm) <- rownames(geo)
dm <- dm / 1000 # kmの単位にする。
dm[1:5, 1:5]
##          Sapporo   Aomori   Morioka   Sendai     Akita
## Sapporo   0.0000 253.8096 373.58228 534.0140 385.85116
## Aomori  253.8096   0.0000 129.30802 283.9589 134.22907
## Morioka 373.5823 129.3080   0.00000 161.1198  90.05514
## Sendai  534.0140 283.9589 161.11979   0.0000 174.19775
## Akita   385.8512 134.2291  90.05514 174.1977   0.00000

距離行列をもとに多次元尺度法を当てはめる。

mds <- cmdscale(dm)
head(mds)
##               [,1]          [,2]
## Sapporo  -887.5166 -381.26148655
## Aomori   -699.6291 -210.60184148
## Morioka  -653.8852  -89.52139326
## Sendai   -539.6094   24.15571884
## Akita    -582.6652 -144.76377589
## Yamagata -502.0629    0.02311809

得点が計算されている。

プロットする。

plot_ly(x = mds[,1], y = mds[,2], text = rownames(mds),
        type = "scatter", mode = "markers")

得点間の関係を確認する。

plot_ly(x = pcs[,1], y = mds[,1], text = rownames(pcs), 
          type = "scatter", mode = "markers")
plot_ly(x = pcs[,2], y = mds[,2], text = rownames(pcs), 
          type = "scatter", mode = "markers")

両者はほぼ直線関係にあることがわかる。両者が完全に直線に乗らない理由は、緯度・経度が完全な直交座標ではないためである。

仮に、緯度経度が直交していて、都市間の距離が緯度・経度座標上のユークリッド距離で表されるとすると（実際はそうではない）、以下のように計算される。

都市間の距離を、緯度経度のユークリッド距離で計算する。

dm.euc <- dist(geo)
mds.euc <- cmdscale(dm.euc)

プロットする。

plot_ly(x = mds.euc[,1], y = mds.euc[,2], text = rownames(mds),
        type = "scatter", mode = "markers")

得点間の関係を確認する。

plot_ly(x = pcs[,1], y = mds.euc[,1], text = rownames(pcs), 
          type = "scatter", mode = "markers")
plot_ly(x = pcs[,2], y = mds.euc[,2], text = rownames(pcs), 
          type = "scatter", mode = "markers")

完全に直線関係になる。主成分分析を行った座標系上で計算されたユークリッド距離をもとに多次元尺度法を適用すると、完全に一致する。

なお、緯度・経度のユークリッド距離として計算されたものと、地理的距離として計算されたものを比較する。

plot(x = as.vector(dm), y = as.vector(as.matrix(dm.euc)),
        xlab = "Euclid-distance based on Latitude/Logitude", ylab = "Geographic distance")

完全には一致しない。これが、主成分分析と多次元尺度法の結果が完全に一致しなかった理由である。