公示地価は不完全ではあるもののパネルデータとして整備することができるため，政策評価などに利用することができるかもしれません．ところが，国土数値情報の公式サイトからzipでDLして解凍し，ArcGISでshapeファイルを読み込んで座標系を変換して，空間結合して最寄り施設までの距離などを計算し，その結果を出力して計量ソフトで読み込み，パネル分析ができるように整形して，…というのはどうにも面倒です．ということでデータのDLから分析までRのみでやってみましょう．

ご注意：このページは RStudio + R Markdown + knit を利用して作成しているため（作成者・保守者には恐縮ですが，黒田のような素人にはどうにも使いづらいのです），パネルデータ作成の例示と分析の足掛かりとなるような必要最小限のコードしか載せていません．前処理が不十分だとか，コマンドの説明が足りないとか，図が汚いとか，地図の凡例が不適切だとか，計量モデルが不適切であるとか，そのような指摘は求めておりません．dplyrで前処理をするべきとか，apply一族と交流せよとか，可視化はggplot2以外認めないとか，その類の宗教的な説法もご遠慮しております．作成にあたって id:yutannihilation さんのブログ記事「RでGISをやるときにはsfパッケージ、という世の中になるらしい。」（2017年1月6日）などを参考にしております．

2017年1月の公開後，分析の追加・微調整を行っております（最終更新：2017年2月）．

Setting

（※ 見出し部分が半端な英語なのは，目次からのページ内ジャンプ機能周りのバグ回避のためです．）

rm(list = ls())
setwd(dir = "c://hoge")
library(curl)
library(sf)

分析期間の設定．

years_int <- c(2010:2016)
years <- substr(x = as.character(years_int), start = 3, stop = 4)  # c("10", "11", ...)

対象範囲を都道府県単位で設定．対応するコードは「都道府県コード」 http://nlftp.mlit.go.jp/ksj/gml/codelist/PrefCd.html などを参照して指定．

prefs <- c("08", "11")  # 茨城県と埼玉県

Download

公示地価のデータを国土数値情報 http://nlftp.mlit.go.jp/ksj/ からDLし，zipを解凍．ちなみに都道府県地価調査のサンプルには地価の単位が異なる「林地」が含まれるので要注意．

f <- tempfile()
for (i in years) {
  for (j in prefs) {
    lprice_path_ij <- paste("http://nlftp.mlit.go.jp/ksj/gml/data/L01/L01-", i, "/L01-", i, "_", j, "_GML.zip", sep = "")
    curl::curl_download(url = lprice_path_ij, destfile = f)
    unzip(zipfile = f, exdir = getwd())
  }
}
unlink(f); rm(f)

リスト型オブジェクトdにデータを代入．年度によってファイルパスなどが異なるので，適宜場合分け．2014年度以前は座標系が指定されていないので，緯度経度であることを指定して読み込む．

d <- list()
for (j in prefs) {
  for (i in 10:11) {
    year_pref <- paste(i, j, sep = "_")
    layer_ij <- paste("L01-", i, "_", j, "-g_LandPrice", sep = "")
    d[[year_pref]] <- sf::st_read(dsn = getwd(), layer = layer_ij, quiet = T, stringsAsFactors = F, crs = "+init=epsg:4612")
  }
  for (i in 12:14) {
    year_pref <- paste(i, j, sep = "_")
    layer_ij <- paste("L01-", i, "_", j, sep = "")
    d[[year_pref]] <- sf::st_read(dsn = getwd(), layer = layer_ij, quiet = T, stringsAsFactors = F, crs = "+init=epsg:4612")
  }
  year15_pref <- paste("15", j, sep = "_")
  d[[year15_pref]] <- sf::st_read(dsn = paste(getwd(), "/L01-15_", j, "_GML", sep = ""), layer = paste("L01-15_", j, sep = ""),  quiet = T, stringsAsFactors = F)
  year16_pref <- paste("16", j, sep = "_")
  d[[year16_pref]] <- sf::st_read(dsn = getwd(), layer = paste("L01-16_", j, sep = ""), quiet = T, stringsAsFactors = F)
}

距離の計算をするために，緯度経度から平面直角座標系（関東はほぼ第9系）に変換する．正式な適用範囲は「平面直角座標系（平成十四年国土交通省告示第九号）」等々を参照．また，価格データが文字列として読み込まれているので，これを数字に変換する．

st_crs(x = d[["16_08"]])  # 確認
## $epsg
## [1] NA
## 
## $proj4string
## [1] "+proj=longlat +ellps=GRS80 +no_defs"
## 
## attr(,"class")
## [1] "crs"
st_is_longlat(x = d[["16_08"]])  # 緯度経度か否か
## [1] TRUE
for (i in years) {
  for (j in prefs) {
    year_pref <- paste(i, j, sep = "_")
    d[[year_pref]] <- st_transform(x = d[[year_pref]], crs = "+init=epsg:2451")  # 座標系の変換
    d[[year_pref]]$L01_006 <- as.numeric(d[[year_pref]]$L01_006)  # 実数の型に変換
  }
}

地価を数値に変換したら，地価分布図が描ける．茨城県ではつくば市が，埼玉県では都心近くの地価が高いことが確認できよう．

d_16 <- rbind.sf(d[["16_08"]], d[["16_11"]])
d_16_lp_f <- cut(x = log(d_16$L01_006), breaks = 10)  # 分割
plot(d_16, col = cm.colors(n = 10)[d_16_lp_f], pch = 19, cex = .5)
legend(x = "topleft", legend = names(table(d_16_lp_f))[10:1], col = cm.colors(n = 10)[10:1], pch = 19, cex = .7)  # 凡例

文字化けしている列はiconv関数で文字コードを変換できる．ここでは市区町村名称を変換してみる．

for (i in c(10:11, 14:16)) {  # 2012年度・2013年度は文字コードが違う模様
  for (j in prefs) {
    year_pref <- paste(i, j, sep = "_")
    d[[year_pref]][, "L01_018"][, 1] <- iconv(d[[year_pref]][, "L01_018"][, 1], from = "utf-8", to = "cp932")
  }
}
head(d[["16_08"]][, "L01_018"][, 1])  # 確認
## [1] 古河 古河 古河 古河 古河 古河
## 44 Levels: かすみがうら つくば つくばみらい ひたちなか 阿見 稲敷 ... 龍ケ崎

Distance to nearest ICs

国土数値情報「高速道路時系列データ」を用いて最寄りIC距離を計算し，道路インフラ整備の影響を検証したいとする．2015年12月31日時点のデータを用いるため，2016年1月1日時点の価格を示している（ことになっている）2016年公示地価まで対応しているものとみなす．実はICに加えてJCTも含まれるので，目的次第ではIC・SIC（スマートIC）に限定したほうがよい．※ 座標系はJGD2011（なので面倒そうな予感が一瞬漂うの）だが，こちらはprjファイルも用意されているのでちょっとラッキー．

DL．

f <- tempfile()
highway_path <- "http://nlftp.mlit.go.jp/ksj/gml/data/N06/N06-15/N06-15_GML.zip"
curl::curl_download(url = highway_path, destfile = f)
unzip(zipfile = f, exdir = getwd())
unlink(f); rm(f)
highway_joint <- sf::st_read(dsn = getwd(), layer = "N06-15_Joint", stringsAsFactors = F, quiet = T)  # ICのpointデータ
iconv(x = head(highway_joint$N06_018), from = "utf-8", to = "cp932")  # IC名
## [1] "南相馬鹿島SASIC" "高岡砺波SIC"     "南砺SIC"         "府中SIC"        
## [5] "豊前"            "上毛PASIC"

座標系を変換．

highway_joint <- st_transform(x = highway_joint, crs = "+init=epsg:2451")

各地点について，全国のICからの最短距離を「最寄りICまでの距離」として付与．ユークリッド距離で単位はメートル．道路ネットワーク距離を利用したい場合はおとなしくArcGIS Network Analystにお世話になるのがよさそう．ひとまずICが供用開始した時点の効果を測るものとする．アナウンスメント効果を考慮するには計画発表のタイミングを使用する必要がある．DIDのフレームで議論する場合にはplacebo test用のplacebo距離も必要となるが，今回は明らかにcommon trendとは考えられないので，以下の分析（FE）では用いない（ちなみにラグを動かしながら実際に計算してみると，地価上昇は供用開始前に収束している様子が観察できる）．

for (i in years_int) {
  i_char <- substr(x = as.character(i), start = 3, stop = 4)
  highway_joint_i <- highway_joint[highway_joint$N06_012 <= i & highway_joint$N06_014 >= i, ]
  highway_joint_i_p <- highway_joint[highway_joint$N06_012 <= (i+2) & highway_joint$N06_014 >= (i+2), ]  # placebo
  for (j in prefs) {
    year_pref <- paste(i_char, j, sep = "_")
    dist_temp <- st_distance(x = d[[year_pref]], y = highway_joint_i)
    d[[year_pref]]$ICdist <- apply(X = dist_temp, MARGIN = 1, FUN = min)
    dist_p_temp <- st_distance(x = d[[year_pref]], y = highway_joint_i_p)  # placebo
    d[[year_pref]]$ICdist_p <- apply(X = dist_p_temp, MARGIN = 1, FUN = min)
  }
}
# hist(d[["16_08"]]$ICdist)  # 確認

IC距離の分布図．

d_16 <- rbind.sf(d[["16_08"]], d[["16_11"]])
d_16_ICdist_f <- cut(x = log(d_16$ICdist), breaks = 10)
plot(d_16, col = cm.colors(n = 10)[d_16_ICdist_f], pch = 19, cex = .5)
legend(x = "topleft", legend = names(table(d_16_ICdist_f))[10:1], col = cm.colors(n = 10)[10:1], pch = 19, cex = .7)
plot(highway_joint, pch = 19, cex = .5, col = 8, add = T)

Cutting sample into sub-sample

サンプルを常磐自動車道周辺に限定する場合．

highway_section <- sf::st_read(dsn = getwd(), layer = "N06-15_HighwaySection", stringsAsFactors = F, quiet = T)  # 路線のlineデータ
highway_section <- st_transform(x = highway_section, crs = "+init=epsg:2451")
# plot(highway_section)
# iconv(x = head(highway_section$N06_007), from = "utf-8", to = "cp932")  # 路線名
joban_index <- iconv(x = highway_section$N06_007, from = "utf-8", to = "cp932") == "常磐自動車道"
# str(highway_section[joban_index, ])  # 確認
highway_joint_joban_is <- st_intersects(x = highway_joint, y = highway_section[joban_index, ], sparse = FALSE)
highway_joint_joban <- highway_joint[apply(X = highway_joint_joban_is, MARGIN = 1, FUN = sum) >= 1, ]
# plot(highway_joint_joban)  # 常磐道のICのみ抽出

2016年度・茨城県を例に，標準地の分布を確認．

d_16_08_joban_dist <- st_distance(x = d[["16_08"]], y = highway_joint_joban)  # 標準地から常磐道ICまでの距離行列
d_16_08_joban <- d[["16_08"]][apply(X = d_16_08_joban_dist, MARGIN = 1, FUN = min) <= 10000, ]  # ICの10km以内の標準地のみ抽出
plot(d_16_08_joban, pch = 19, cex = .5, col = 8)
plot(highway_section[joban_index, ], add = T, pch = 19)  # 路線
plot(highway_joint_joban, add = T, pch = 19, col = 2)  # IC点

Shaping panel data

座標を地点IDとして利用．

for (i in years) {
  for (j in prefs) {
    year_pref <- paste(i, j, sep = "_")
    d_coord_ij <- do.call(rbind, unclass(st_geometry(d[[year_pref]])))
    # d[[year_pref]]$coordx <- d_coord_ij[,1]
    # d[[year_pref]]$coordy <- d_coord_ij[,2]
    d[[year_pref]]$coord <- paste(d_coord_ij[,1], d_coord_ij[,2], sep = "_")
  }
}

まとめる．

pdata_cols <- c("coord", "L01_005", "L01_001", "L01_002", "L01_006", "L01_018", "ICdist", "ICdist_p")
d_merged <- NULL
for (i in years) {
  for (j in prefs) {
    year_pref <- paste(i, j, sep = "_")
    y_i <- d[[year_pref]][, pdata_cols]
    d_merged <- rbind(d_merged, y_i)
  }
}

str(d_merged)  # 確認
## 'data.frame':    13751 obs. of  8 variables:
##  $ coord   : chr  "-11767.2626678934_21472.1520597889" "-11748.2879866905_20736.1837188635" "-11703.8787962622_22328.9479794719" "-11617.7575261593_21633.1646811712" ...
##  $ L01_005 : chr  "2010" "2010" "2010" "2010" ...
##  $ L01_001 : chr  "000" "000" "000" "005" ...
##  $ L01_002 : chr  "003" "004" "002" "004" ...
##  $ L01_006 : num  56600 49700 48600 71400 84200 54800 76000 40200 48800 30500 ...
##  $ L01_018 : chr  "古河" "古河" "古河" "古河" ...
##  $ ICdist  : num  10011 9494 10700 10231 10350 ...
##  $ ICdist_p: num  10011 9494 10700 10231 10350 ...
head(table(d_merged$L01_018))  # 自治体ごとの地点数
## 
## かすみがうら   さいたま西 さいたま中央   さいたま南   さいたま北 
##           82           62          124          148          190 
##   さいたま緑 
##          136
table(table(d_merged$coord))  # 重複度ごとの地点数
## 
##    1    2    3    4    5    6    7 
##  392   96  566  620   59   42 1206
sum(table(table(d_merged$coord)))  # 地点数 length(unique(d_merged$coord))
## [1] 2981

Fixed-effects model

パッケージ plm を用いて固定効果分析を行う．Key variableはIC距離．両対数型（log-log）で定式化すれば弾力性を推定することになる（はず．一応は）． \[ \log (Land Price)_{i, t} = \alpha_{i} + \alpha_{t} + \beta * \log (ICdist)_{i, t} + \epsilon_{i, t} \]

fm00 <- log(L01_006) ~ log(ICdist)

推定に時間がかかるのでサンプルの一部を切り出す．

d_merged_sub <- d_merged[d_merged$L01_005 %in% 2014:2016, ]
str(d_merged_sub)
## 'data.frame':    5570 obs. of  8 variables:
##  $ coord   : chr  "-11767.2626678934_21472.1520597889" "-11706.5925763133_20708.1678779883" "-11703.8787962622_22328.9479794719" "-11461.5750079171_21651.8754091764" ...
##  $ L01_005 : chr  "2014" "2014" "2014" "2014" ...
##  $ L01_001 : chr  "000" "000" "000" "005" ...
##  $ L01_002 : chr  "003" "004" "002" "001" ...
##  $ L01_006 : num  52800 44400 43900 75300 71500 35800 43900 29700 36000 70000 ...
##  $ L01_018 : chr  "古河" "古河" "古河" "古河" ...
##  $ ICdist  : num  10011 9504 10700 10350 10723 ...
##  $ ICdist_p: num  10011 9504 10700 10350 10723 ...

どこでハンドリングを誤ったのか地点と時点のペアに重複が存在するので，繰り返しが出現した部分は除去する．気分が悪いが，とりあえず姑息なエラー回避で凌ぐ．

d_coord_year <- paste(d_merged_sub$coord, d_merged_sub$L01_005)
d_merged_sub <- d_merged_sub[!duplicated(d_coord_year), ]
str(d_merged_sub)
## 'data.frame':    5569 obs. of  8 variables:
##  $ coord   : chr  "-11767.2626678934_21472.1520597889" "-11706.5925763133_20708.1678779883" "-11703.8787962622_22328.9479794719" "-11461.5750079171_21651.8754091764" ...
##  $ L01_005 : chr  "2014" "2014" "2014" "2014" ...
##  $ L01_001 : chr  "000" "000" "000" "005" ...
##  $ L01_002 : chr  "003" "004" "002" "001" ...
##  $ L01_006 : num  52800 44400 43900 75300 71500 35800 43900 29700 36000 70000 ...
##  $ L01_018 : chr  "古河" "古河" "古河" "古河" ...
##  $ ICdist  : num  10011 9504 10700 10350 10723 ...
##  $ ICdist_p: num  10011 9504 10700 10350 10723 ...

plm関数では引数indexによる指定を省略すると，データフレームの1列目をindividualとして，2列目をtimeとしてそれぞれ認識する．

library(plm, quietly = T)
plm00 <- plm(formula = fm00, data = d_merged_sub, effect = "twoways", model = "within")
summary(plm00)
## Twoways effects Within Model
## 
## Call:
## plm(formula = fm00, data = d_merged_sub, effect = "twoways", 
##     model = "within")
## 
## Unbalanced Panel: n=2030, T=1-3, N=5569
## 
## Residuals :
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -0.074300 -0.004970  0.000168  0.004300  0.077100 
## 
## Coefficients :
##               Estimate Std. Error t-value Pr(>|t|)
## log(ICdist) -0.0014773  0.0022549 -0.6551   0.5124
## 
## Total Sum of Squares:    0.72528
## Residual Sum of Squares: 0.72519
## R-Squared:      0.00012136
## Adj. R-Squared: -0.57447
## F-statistic: 0.429185 on 1 and 3536 DF, p-value: 0.51243
# summary(plm::fixef(plm00, type = "dmean"))  # individual FE
summary(plm::fixef(plm00, type = "dmean", effect = "time"))
##        Estimate Std. Error t-value Pr(>|t|)
## 2014 -0.0010482  0.0185290 -0.0566   0.9549
## 2015 -0.0044926  0.0184473 -0.2435   0.8076
## 2016  0.0055408  0.0184354  0.3006   0.7638

サンプルが2012-13年を含む場合は，その部分だけ「西暦 201x」のようになっているが，気にしない．

ちなみに，1時点分しかない地点を除外するには以下のようにサンプルを切り出す．

d_merged_sub_dup <- names(table(d_merged_sub$coord))[table(d_merged_sub$coord) >= 2]
d_merged_sub2 <- d_merged_sub[d_merged_sub$coord %in% d_merged_sub_dup, ]
plm01 <- plm(formula = fm00, data = d_merged_sub, effect = "twoways", model = "within")
summary(plm01)
## Twoways effects Within Model
## 
## Call:
## plm(formula = fm00, data = d_merged_sub, effect = "twoways", 
##     model = "within")
## 
## Unbalanced Panel: n=2030, T=1-3, N=5569
## 
## Residuals :
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -0.074300 -0.004970  0.000168  0.004300  0.077100 
## 
## Coefficients :
##               Estimate Std. Error t-value Pr(>|t|)
## log(ICdist) -0.0014773  0.0022549 -0.6551   0.5124
## 
## Total Sum of Squares:    0.72528
## Residual Sum of Squares: 0.72519
## R-Squared:      0.00012136
## Adj. R-Squared: -0.57447
## F-statistic: 0.429185 on 1 and 3536 DF, p-value: 0.51243
coef(plm00)[1] == coef(plm01)[1]  # 一致するはず
## log(ICdist) 
##        TRUE

plm::plmを使わなくてもlm関数で推定可能．

fm01 <- update.formula(old = fm00, new = . ~ . + coord + L01_005)
lm01 <- lm(formula = fm01, data = d_merged_sub)
head(summary(lm01)$coefficients, 2)
##                 Estimate  Std. Error     t value Pr(>|t|)
## (Intercept) 11.044533724 0.019721363 560.0289312 0.000000
## log(ICdist) -0.001477262 0.002254943  -0.6551218 0.512432

上記は誤差項に対して分散均一を仮定して標準誤差を計算しているため，この妥当性を確認する．p値が大きい時（たとえば p > 0.05）は「帰無仮説：均一分散」を棄却できないが，p値が小さい場合は不均一分散が示唆される．

library(lmtest, quietly = T)
lmtest::bptest(plm00)  # Breusch-Pagan test for heteroskedasticity
## 
##  studentized Breusch-Pagan test
## 
## data:  plm00
## BP = 27.876, df = 1, p-value = 1.293e-07

地点ごとのばらつきの違いを反映した標準誤差に基づいてp値を計算するために，clustered stanrard errorを用いる．コードはAjay Shah先生によるブログ記事「Sophisticated clustered standard errors using recent R tools」（2016年6月15日）を参考にしている．

lmtest::coeftest(x = plm00, vcov = plm::vcovHC(x = plm00, type = "sss", cluster = "group"))
## 
## t test of coefficients:
## 
##               Estimate Std. Error t value Pr(>|t|)
## log(ICdist) -0.0014773  0.0012679 -1.1652    0.244

残差の空間分布を確認（七面倒なハンドリングをしているが，なんとかならないものだろうか）．

plm00_2016 <- attr(plm00$model, "index")$L01_005 == "2016"  # logical: 2016年かどうか
plm00_coord <- do.call(rbind, strsplit(x = as.character(attr(plm00$model, "index")$coord), split = "_"))
plm00_x <- as.numeric(x = plm00_coord[, 1])  # 全期間
plm00_y <- as.numeric(x = plm00_coord[, 2])
plm00_x_2016 <- plm00_x[plm00_2016]  # 2016年分
plm00_y_2016 <- plm00_y[plm00_2016]

plm00_resid_16_f <- cut(x = residuals(plm00)[plm00_2016], breaks = 10)
plot(x = plm00_x_2016, y = plm00_y_2016, col = cm.colors(n = 10)[plm00_resid_16_f], axes = F, ann = F, pch = 19, cex = .5)
legend(x = "topleft", legend = names(table(plm00_resid_16_f))[10:1], col = cm.colors(n = 10)[10:1], pch = 19, cex = .7)

# plot(highway_joint, add = T, pch = 19, cex = .5, col = 8)  # IC等

Common trend assumption

の予定だったが，今回のケースではcommon trendのはずがないので，以下はplacebo testのイメージということで． 

plot(x = d_merged_sub$ICdist, y = d_merged_sub$ICdist_p, pch = 20, cex = .3)  # treatment vs. placebo treatment

table(d_merged_sub$ICdist == d_merged_sub$ICdist_p)
## 
## FALSE  TRUE 
##   140  5429

\[ \log (Land Price)_{i, t} = \alpha_{i} + \alpha_{t} + \beta * \log (ICdist)_{i, t+2} + \epsilon_{i, t} \]

fm01_p <- log(L01_006) ~ log(ICdist_p) + coord + L01_005
lm01_p <- lm(formula = fm01_p, data = d_merged_sub)
head(summary(lm01_p)$coefficients, 2)
##                 Estimate  Std. Error   t value Pr(>|t|)
## (Intercept)   13.9666287 0.047881605 291.69091        0
## log(ICdist_p) -0.3909362 0.005553463 -70.39502        0
# summary(lm01_p)$coefficients["log(ICdist_p)", ]
coef(lm01)[2] == coef(lm01_p)[2]  # 異なるはず
## log(ICdist) 
##       FALSE

Cross-term

固定効果モデルでは残差に空間パターンが残るため，地域 \(r = r(i)\) （市区町村）と時点の交差項を含める． \[ \log (Land Price)_{i, t} = \alpha_{i} + \alpha_{t} + \alpha_{r, t} + \beta * \log (ICdist)_{i, t} + \epsilon_{i, t} \]

fm03 <- update.formula(old = fm01, new = . ~ . + coord + L01_005 + L01_018:L01_005)
lm03 <- lm(formula = fm03, data = d_merged_sub)
head(summary(lm03)$coefficients, 2)
##                 Estimate  Std. Error     t value  Pr(>|t|)
## (Intercept) 1.104167e+01 0.020704299 533.3031411 0.0000000
## log(ICdist) 4.628526e-04 0.002597004   0.1782256 0.8585567

ナイーブにclustered SEを計算．

library(multiwayvcov)
lm03_cse_i <- coeftest(x = lm03, vcov = cluster.vcov(model = lm03, cluster = d_merged_sub$coord))
lm03_cse_i[1:2, ]  # clustered SE: individual clustering
##                 Estimate  Std. Error    t value  Pr(>|t|)
## (Intercept) 1.104167e+01 0.023856970 462.827747 0.0000000
## log(ICdist) 4.628526e-04 0.002937385   0.157573 0.8748029
lm03_cse_m <- coeftest(x = lm03, vcov = cluster.vcov(model = lm03, cluster = d_merged_sub$L01_018))
lm03_cse_m[1:2, ]  # clustered SE: municipality clustering
##                 Estimate  Std. Error      t value  Pr(>|t|)
## (Intercept) 1.104167e+01 0.042062452 262.50651720 0.0000000
## log(ICdist) 4.628526e-04 0.005512457   0.08396485 0.9330895

残差の空間分布を確認．

lm03_16 <- lm03$model$L01_005 == "2016"  # logical: 2016年かどうか
lm03_coord <- do.call(rbind, strsplit(x = lm03$model$coord, split = "_"))
lm03_x <- as.numeric(x = lm03_coord[, 1])  # 全期間
lm03_y <- as.numeric(x = lm03_coord[, 2])
lm03_x_2016 <- lm03_x[lm03_16]  # 2016年分
lm03_y_2016 <- lm03_y[lm03_16]

lm03_resid_16_f <- cut(x = residuals(lm03)[lm03_16], breaks = 10)
plot(x = lm03_x_2016, y = lm03_y_2016, col = cm.colors(n = 10)[lm03_resid_16_f], pch = 19, cex = .5)
legend(x = "topleft", legend = names(table(lm03_resid_16_f))[10:1], col = cm.colors(n = 10)[10:1], pch = 19, cex = .7)

（以上）