よくデータ解析の例で使われるデータに boston housing data があるが， 似たようなことを我らが豊中でやってみようと思い，賃貸物件サイトから豊中市の物件のデータを取ってみた． 今回はこのデータについて，簡単に説明する．

データについて

データとして，2015年2月14日時点でのいい部屋ネットの「豊中市」でヒットした物件すべてのデータを取得した． 実際に使うデータは，価格・位置が重複する物件（同一建物で部屋違い）を省いたもので，省いた後のサンプル数は 6,863．

変数は 153 種ある．

• 住所 / 最寄り駅 / 駅からの所要時間 / バス利用かどうか
• 価格（月額）/ 管理費 / 敷金 / 礼金 / 償却金
• 間取り / 占有面積 / 所在階 / 方角
• 建物の総階数 / 戸数 / 築年数 / 完成年月 / 建物種別 / 構造
• 緯度 / 経度
• 駐車場の有無 / 敷地内駐車場かどうか / 駐車場価格
• 取扱店舗 / 情報公開日 / 契約種別 / 契約期間
• その他，設備（オートロックなど）の有無を表す二値変数 116 種

業者によって入力されている項目はまちまちなので，データにはそこそこの数の欠損が含まれる． とくに最後の設備の有無については，実際には設備があっても入力されていなければ 0 なので，注意が必要．

地図上での賃貸価格分布

まずは地図上で賃貸価格の分布を見てみる． RgoogleMaps パッケージを使うと，こんな感じで簡単に地図のデータプロットを作成できる．

library(RgoogleMaps)
center <- c(mean(range(dat$lat)), mean(range(dat$lon)))
map <- GetMap(center, zoom=13)
PlotOnStaticMap(map, lat=dat$lat, lon=dat$lon, pch=3, col=col)

で，実際のプロットがこちら．色は価格の対数でスケールさせている．

この図から，豊中市の南北の差が極めてはっきり読み取れる． 千里中央・緑丘（右上のあたり）は高級住宅地だけあって高額物件が集中していて，安い物件はまったく見当たらない． 中央あたりの豊中駅周辺・岡町まで来ると高額物件と安い物件が混ざった状態になり， 曽根を過ぎるあたりから高額物件の青い点が減っている． 蛍池・伊丹空港あたりに安い物件が集中しているのは，近くに阪大豊中キャンパスがあるせいだろうか．

とりあえず地図上にプロットしてみただけだが，地域ごとの偏りがはっきり見えていて，なかなか面白い．

価格分布のヒストグラム

次に，価格のヒストグラムをとってみる．

breaks <- c(1.25, 2.5, 5, 10, 20, 40) * 10000
g <- ggplot(dat, aes(x=price))
g <- g + scale_x_continuous(trans="log", breaks=breaks)
g <- g + geom_histogram(aes(y=..density..)) 
g <- g + stat_function(fun=function(x) dnorm(x, mean(log(dat$price)), sd(log(dat$price))))
g

横軸は対数軸で，青い線は標本平均と不偏分散から作った正規分布の密度関数． これを見ると，多少頭が抜け出ているものの，全体として非常にきれいに対称に分布していることがわかる． この図を見る限りでは，賃貸価格は対数正規分布に従うとみてよさそうである．

建築年のヒストグラム

あとは面白そうなところで，建築年のヒストグラムを見てみる．

自然に考えると古い方から新しい方へだんだん増加するように思えるが，実際には山や谷が強く出ている． 物件の多い時期を見ていくと，1990年前後はバブル期，1996年は阪神大震災の影響が思い当たる． 物件の少ない方では，1981年と2000年に建築基準法改正があったが，なぜ改正されると物件が減るかというのは説明できないので，本当にそれが原因かはよく分からない． なお，このデータは「売り出し中の」物件データなので，物件全体と比べて多少のバイアスが入っていることに注意．

And more…

今回はデータの概要を紹介したが，次回からはこのデータを使っていろいろな分析を試していきたいと思う． 賃貸価格予測のほか，部屋選びで重視されるポイントを調べたり，いわゆる「いわくつき物件」の判定ができたら面白そうだ．

今回使用したスクリプト

library(ggplot2) library(RgoogleMaps) options(scipen=1)

dat.column <- read.csv("column.csv", head=FALSE, fileEncoding="UTF-8", colClasses="character") dat <- read.csv("data.csv", fileEncoding="UTF-8", colClasses=dat.column[, 3]) dat <- dat[!duplicated(cbind(dat$price, dat$lat, dat$lon)), ] dat$price[dat$id == "300056150013621013622"] <- NA # 価格の桁間違い? dat$area[dat$id %in% c("300054130041305041305", "300105420003450004117")] <- NA # 面積の桁間違い? dat$latitude[dat$id == "300059820006005006005"] <- NA # 緯度不正 dat$longitude[dat$id == "300059820006005006005"] <- NA # 経度不正 dim(dat) ## ---- ## pricehist mean.log.price <- mean(log(dat$price), na.rm=TRUE) sd.log.price <- sd(log(dat$price), na.rm=TRUE) breaks <- c(1.25, 2.5, 5, 10, 20, 40) * 10000 g <- ggplot(dat, aes(x=price)) g <- g + scale_x_continuous( trans="log", breaks=breaks, labels=breaks / 10000, name="価格 (万円)") g <- g + theme(panel.grid.minor=element_blank()) g <- g + geom_histogram(aes(y=..density..), col=hsv(0.6, 1, 0.75), fill=hsv(0.6, 0.5, 1)) g <- g + stat_function(fun=function(x) dnorm(x, mean.log.price, sd.log.price), size=1, col="blue") g ## ---- ## builthist dat.building <- dat[!duplicated(cbind(dat$lat, dat$lon)), ] year <- as.integer(substr(dat.building$built, 1, 4)) g <- ggplot(data.frame(year=na.omit(year)), aes(x=year)) g <- g + scale_x_continuous(breaks=seq(1950, 2010, by=10), name="年") g <- g + geom_histogram(col=hsv(0.6, 1, 0.75), fill=hsv(0.6, 0.5, 1), binwidth=1) g ## ---- ## map dat.map <- na.omit(dat[, c("price", "latitude", "longitude")]) center <- c(mean(range(dat.map$lat)), mean(range(dat.map$lon))) map <- GetMap(center, zoom=13) lv2col <- function(lv, ...) hsv(pmin(pmax(lv, 0), 1) * 0.75, 1, 1, ...) q <- quantile(log(dat.map$price), c(0.05, 0.95), na.rm=TRUE) price.level <- (log(dat.map$price) - q[1]) / (q[2] - q[1]) PlotOnStaticMap(map, lat=dat.map$lat, lon=dat.map$lon, pch=3, col=lv2col(price.level)) legend.level <- seq(0, 1, by=0.2) legend("bottomleft", legend=round(exp(legend.level * (q[2] - q[1]) + q[1])), fill=lv2col(legend.level)) [/code]