前回の記事では，豊中市の賃貸物件データの紹介をした． 今回の記事では，どの変数が賃貸価格に影響を及ぼしているかを調べるために，賃貸データに対して回帰分析を適用してみる． 比較のために，OLS，Lasso，Adaptive Lasso の3つの手法を実行して結果を見ることにした．
使用するデータ
目的変数には対数をとった賃貸価格， 説明変数には文字列・カテゴリカル変数を除いたすべての変数と，一部のカテゴリカル変数をダミー変数化して使う． 変数はあらかじめすべて標準化し，NA を含むサンプルは除外している．
分析にかける前に，相関の高い変数がないかチェックしておく．
cor.dat <- cor(dat.explanner) large.cor <- matrix(abs(cor.dat) > 0.7, nrow(cor.dat)) & upper.tri(cor.dat) large.cor.index <- arrayInd(which(large.cor), .dim=dim(large.cor)) dat.large.cor <- data.frame( cor=cor.dat[large.cor], col=matrix(colnames(dat.explanner)[large.cor.index], nrow(large.cor.index)), text=apply(large.cor.index, 2, function(r) varnames[colnames(dat.explanner)[r]]) ) dat.large.cor[order(abs(dat.large.cor$cor), decreasing=TRUE), ]
cor col.1 col.2 text.1 text.2
9 0.8805938 E7 E89 冷房 暖房機
1 0.8355318 num.room area 部屋数 面積
5 0.8353020 E48 E72 室内洗濯機置き場 シャワー
11 0.7984594 E84 E92 敷地内ごみ置場 専用ゴミ置場
10 -0.7723538 E81 E92 給湯 専用ゴミ置場
8 0.7694277 E53 E88 電話機 テレビ
7 0.7634514 E6 E87 可動間仕切り アプローチライト
6 -0.7285695 E72 E81 シャワー 給湯
3 0.7254785 type.mansion struct.concrete マンション コンクリート建築
2 0.7152430 floor num.floor 所在階 総階数
4 0.7013576 E21 E72 バルコニー シャワー
気づいたこと:
- 冷房と暖房は両方ついた部屋が多いとか，（アパートや一戸建て以外という意味での）マンションはコンクリート建築が多いとかはなるほどなあと思った．
- 可動間仕切りとアプローチライトの相関の高さは何だ？どちらもサンプル数少ないから影響少なそうだけど……． 可動間仕切りって障子とかですかね．
- シャワーがやたら出てくる．給湯とかバルコニーとか，一定以上の水準ならほとんどあるものが拾われている？
上の結果を踏まえて，次のように変数を変更した．
- 冷房と暖房機は OR をとって新たな変数「冷暖房」にまとめる．
- 敷地内ごみ置き場は専用ゴミ置き場にカウント．
- 部屋数，シャワー，電話機，可動間仕切り，マンション，総階数は除外．
この変更の結果，相関係数が 0.7 以上の変数はなくなった． 本当は 0.5 くらいまではチェックしたほうがよいのだろうが，とりあえずこのくらいで．
OLS による賃貸価格の回帰分析
まずはオーソドックスに OLS（最小二乗法）による線形回帰から． ダミー変数の係数比較に意味があるかとか細かい議論は置いておいて，とりあえずざっくり lm にかけて係数を比較してみる． 変数が多いので，BIC による変数減少法を使って変数選択を行う．
biclm.dat <- step(lm(log.price ~ ., data=dat.lm), k=log(nrow(dat.lm))) summary(biclm.dat)
Call:
lm(formula = log.price ~ time + L + D + area + age + floor +
type.house + struct.mokuzo + struct.concrete + C11 + C15 +
C16 + E2 + E5 + E10 + E12 + E13 + E17 + E20 + E25 + E26 +
E28 + E30 + E34 + E37 + E41 + E45 + E46 + E47 + E48 + E55 +
E56 + E57 + E60 + E67 + E81 + E90 + E96 + E7E89, data = dat.lm)
Residuals:
Min 1Q Median 3Q Max
-2.3516 -0.2039 -0.0011 0.2066 1.9305
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.717e-15 4.427e-03 0.000 1.000000
time -1.030e-01 4.749e-03 -21.687 < 2e-16 ***
L 6.571e-02 6.956e-03 9.446 < 2e-16 ***
D 7.672e-02 6.456e-03 11.883 < 2e-16 ***
area 6.505e-01 7.840e-03 82.968 < 2e-16 ***
age -2.129e-01 6.568e-03 -32.421 < 2e-16 ***
floor 4.044e-02 5.056e-03 7.998 1.50e-15 ***
type.house 5.932e-02 5.983e-03 9.914 < 2e-16 ***
struct.mokuzo -3.029e-02 5.691e-03 -5.322 1.06e-07 ***
struct.concrete 3.153e-02 6.142e-03 5.133 2.93e-07 ***
C11 -1.471e-02 4.490e-03 -3.275 0.001061 **
C15 1.792e-02 4.590e-03 3.903 9.60e-05 ***
C16 -2.026e-02 4.516e-03 -4.488 7.34e-06 ***
E2 4.342e-02 4.915e-03 8.834 < 2e-16 ***
E5 1.025e-01 5.127e-03 19.992 < 2e-16 ***
E10 -2.414e-02 4.873e-03 -4.954 7.48e-07 ***
E12 -1.554e-02 4.716e-03 -3.296 0.000986 ***
E13 3.405e-02 5.224e-03 6.518 7.72e-11 ***
E17 3.349e-02 5.261e-03 6.365 2.09e-10 ***
E20 -1.462e-02 4.871e-03 -3.002 0.002694 **
E25 2.914e-02 6.119e-03 4.763 1.96e-06 ***
E26 1.595e-02 4.661e-03 3.422 0.000626 ***
E28 -2.136e-02 5.694e-03 -3.752 0.000177 ***
E30 3.262e-02 5.399e-03 6.042 1.61e-09 ***
E34 3.314e-02 5.933e-03 5.585 2.44e-08 ***
E37 3.021e-02 5.709e-03 5.293 1.25e-07 ***
E41 3.301e-02 5.469e-03 6.035 1.68e-09 ***
E45 1.478e-02 4.694e-03 3.148 0.001652 **
E46 2.084e-02 4.675e-03 4.458 8.41e-06 ***
E47 2.099e-02 4.754e-03 4.414 1.03e-05 ***
E48 -4.039e-02 7.818e-03 -5.166 2.47e-07 ***
E55 1.792e-02 5.628e-03 3.185 0.001456 **
E56 4.111e-02 5.503e-03 7.471 9.09e-14 ***
E57 2.226e-02 5.320e-03 4.185 2.90e-05 ***
E60 1.801e-02 5.611e-03 3.209 0.001337 **
E67 1.705e-02 4.501e-03 3.789 0.000153 ***
E81 -5.480e-02 7.447e-03 -7.359 2.10e-13 ***
E90 1.648e-02 4.678e-03 3.522 0.000431 ***
E96 3.235e-02 5.768e-03 5.609 2.12e-08 ***
E7E89 -1.645e-02 4.931e-03 -3.335 0.000858 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.347 on 6104 degrees of freedom Multiple R-squared: 0.8804, Adjusted R-squared: 0.8796 F-statistic: 1152 on 39 and 6104 DF, p-value: < 2.2e-16
棒グラフのエラーバーは t 検定による 95% 信頼区間である．
まず注目すべきは決定係数の高さ． 説明変数が多いとはいえ，決定係数 0.880 というのは非常によいフィッティングだといえるだろう．
また，変数選択に BIC を使っているので，変数の数は 39 とだいぶ絞られていている．
値を見てみるとすべての変数で 0.3% 以下になっているので，
よく言われている通り BIC はかなり保守的な選択をしていると言える．
具体的に係数の大小を見ていくと，面積，築年数，駅からの所要時間などの影響が強いのはまあそうだろうなという感じ． 全体的に，係数が正の部分は普通のアパートならなさそうな高級感ある要素が並んでいてわかりやすい． 一方で負の部分は価格の高低にかかわらず有無が決まりそうな変数ではあるが，なぜ値段が下がるのかちょっと解釈しづらい感がある．
面白いと思ったのは，バス・トイレ別のセパレートタイプかどうかが所要時間と同じくらい評価されていること． そういえば自分が今の部屋を選んだ時もセパレートは条件だったし，周りに聞いてみるとやはりセパレートがいいとの 意見が大勢だったので，なんとなく感じていたセパレートの重要さがデータによって裏付けられたといえる．
ところで，事務所利用可とか防犯シャッターがあると値段下がるのって，やっぱりアレなんだろうか．
Lasso による賃貸価格の回帰分析
続いて Lasso による回帰を試す． 推定には glmnet パッケージを使い，チューニングパラメータは BIC によって決定する．
X <- as.matrix(dat.explanner)
y <- dat.lm$log.price
glmnet.dat <- glmnet(X, y)
cri <- mapply(function(lambda, df) {
beta <- as.numeric(coef(glmnet.dat, lambda))
resid <- drop(y - cbind(rep(1, nrow(X)), X) %*% beta)
log(var(resid)) + log(nrow(X)) * df / nrow(X)
}, glmnet.dat$lambda, glmnet.dat$df)
best.i <- which.min(cri)
coef.glmnet.dat <- drop(coef(glmnet.dat, glmnet.dat$lambda[best.i]))[-1]
選択された変数の数は 51 と増えているが，上位下位は lm とほとんど変わらない印象． とはいえ，lm で変数減少法をすると 5 分以上かかるのに対して，glmnet では 1 秒で終わってしまうので， 実用性で見れば Lasso での変数選択は変数減少法よりずっと使い勝手がよいといえる．
Adaptive Lasso による賃貸価格の回帰分析
上の結果を助教殿に見せたら「Adaptive Lasso でもやってみてよ」と言われたので， 助教殿謹製の msgps パッケージを使って Adaptive Lasso を試してみることにした．
通常の Lasso では，L1 罰則項つき対数尤度に基づく目的関数
![\[ L_\text{Lasso} = \| \bm{y} - \sum_{j=1}^{p} \bm{x}_{j} \beta_j \|^2 + \lambda \sum_{j=1}^{p} | \beta_j | \]](http://www.hnagata.net/wp/wp-content/ql-cache/quicklatex.com-461a0e6b9045bf21301117c72eede9df_l3.png "Rendered by QuickLaTeX.com")
を最小化するが，Adaptive Lasso ではペナルティ項にウェイトを付けた目的関数
![\[ L_\text{ALasso} = \| \bm{y} - \sum_{j=1}^{p} \bm{x}_{j} \beta_j \|^2 + \lambda \sum_{j=1}^{p} \frac{|\beta_j|}{|\hat\beta_j|^\gamma} \]](http://www.hnagata.net/wp/wp-content/ql-cache/quicklatex.com-767dff99e97e80ed2bc2dbdadcf96c52_l3.png "Rendered by QuickLaTeX.com")
を最小化する．
は 
の一致推定量ならなんでもよく，
たいていは OLS で求められる最尤推定量が使われるので今回もそれに従う．
また，
は2つ目のチューニングパラメータでウェイトの強さを決定するが，
ここでは msgps パッケージのデフォルトにならって 
を使う．
Lasso は変数選択についての一致性を持たないが，Adaptive Lasso は適当な条件の下で変数選択の一致性を持つようになっているという意味で， Adaptive Lasso は Lasso の一種の改良になっている．
とりあえず Adaptive Lasso をやってみたが，今回は lm で十分うまくいっている状況なので， 結局 Lasso も Adaptive Lasso も精度的にはほとんど変わらなかった． サンプル数が少ない状況ではもっと差が出るのだろうか．
今回使用したスクリプト
library(ggplot2) library(glmnet) library(msgps) options(scipen=1)
dat.column <- read.csv("column.csv", head=FALSE, fileEncoding="UTF-8", colClasses="character") varnames <- dat.column[, 2] names(varnames) <- dat.column[, 1] 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) dat.lm <- cbind( log.price = log(dat$price), dat[, c("time", "num.room", "S", "L", "D", "K", "area", "age", "floor", "num.floor", "days.published")], type.mansion = as.numeric(dat$type == "マンション"), type.house = as.numeric(grepl("一戸建て", dat$type)), struct.mokuzo = as.numeric(dat$struct == "木造"), struct.concrete = as.numeric(dat$struct == "ALC" | grepl("コンクリート", dat$struct)), dir.N = as.numeric(dat$dir == "北"), dir.NE = as.numeric(dat$dir == "北東"), dir.E = as.numeric(dat$dir == "東"), dir.SE = as.numeric(dat$dir == "南東"), dir.S = as.numeric(dat$dir == "南"), dir.SW = as.numeric(dat$dir == "南西"), dir.W = as.numeric(dat$dir == "西"), parking = as.numeric(dat$parking != "なし"), dat[, paste0("C", 1:17)], dat[, paste0("E", (1:99)[-39])] ) dat.lm <- dat.lm[, -which(colnames(dat.lm) == "E89")] # 暖房機: E7冷房と高相関 (0.881) dat.lm <- dat.lm[, -which(colnames(dat.lm) == "num.room")] # area と高相関 (0.836) dat.lm <- dat.lm[, -which(colnames(dat.lm) == "E72")] # シャワー: E48洗濯機置き場と高相関 (0.835) dat.lm$E92 <- as.numeric(dat.lm$E92 | dat.lm$E84) dat.lm <- dat.lm[, -which(colnames(dat.lm) == "E84")] # 敷地内ゴミ置き場: E92専用ゴミ置き場と統合 (0.798) dat.lm <- dat.lm[, -which(colnames(dat.lm) == "E53")] # 電話機: E88テレビと高相関 (0.769) dat.lm <- dat.lm[, -which(colnames(dat.lm) == "E6")] # 可動間仕切り: E87アプローチライトと高相関 (0.763) dat.lm <- dat.lm[, -which(colnames(dat.lm) == "type.mansion")] # struct.concrete と高相関 (0.725) dat.lm <- dat.lm[, -which(colnames(dat.lm) == "num.floor")] # floor と高相関 (0.715) dat.lm <- as.data.frame(scale(na.omit(dat.lm))) dim(dat.lm) varnames <- c( varnames, log.price = "log(価格)", type.mansion = "マンション", type.house = "一戸建て", struct.mokuzo = "木造", struct.concrete = "コンクリート建築", dir.N = "北向き", dir.NE = "北東向き", dir.E = "東向き", dir.SE = "南東向き", dir.S = "南向き", dir.SW = "南西向き", dir.W = "西向き" ) ## ---- ## check corelations cor.dat <- cor(dat.lm) large.cor <- matrix(abs(cor.dat) > 0.7, nrow(cor.dat)) & upper.tri(cor.dat) large.cor.index <- arrayInd(which(large.cor), .dim=dim(large.cor)) cbind( matrix(colnames(dat.lm)[large.cor.index], nrow(large.cor.index)), cor.dat[large.cor] ) ## ---- ## plot function my.barplot.lm <- function(var, est, err=NULL) { var <- factor(varnames[var], levels=varnames[var]) ggdat <- data.frame(var=var, est=est, sign=est>=0) if (!is.null(err)) ggdat$err = err g <- ggplot(ggdat, aes(x=var, y=est, colour=sign, fill=sign)) g <- g + coord_flip() g <- g + scale_colour_manual(values=c(hsv(0, 1, 0.75), hsv(0.6, 1, 0.75))) g <- g + scale_fill_manual(values=c(hsv(0, 0.5, 1), hsv(0.6, 0.5, 1))) g <- g + theme(legend.position="none") g <- g + scale_x_discrete(name="") g <- g + scale_y_continuous(name="推定値") g <- g + geom_bar(stat="identity") if (!is.null(ggdat$err)) { g <- g + geom_errorbar(aes(ymin=est-err, ymax=est+err), width=0.3) } g } ## ---- ## lm if (!file.exists("aiclm.Rdata")) { aiclm.dat <- step(lm(log.price ~ ., data=dat.lm)) save(aiclm.dat, file="aiclm.Rdata") } else { load("aiclm.Rdata") } coef.aiclm.dat <- summary(aiclm.dat)$coef[-1, ] coef.aiclm.dat <- coef.aiclm.dat[order(coef.aiclm.dat[, 1]), ] my.barplot.lm( rownames(coef.aiclm.dat), coef.aiclm.dat[, 1], err=coef.aiclm.dat[, 2] * qnorm(0.975)) ggsave("aiclm.svg", width=6, height=8, family="Meiryo") if (!file.exists("biclm.Rdata")) { biclm.dat <- step(lm(log.price ~ ., data=dat.lm), k=log(nrow(dat.lm))) save(biclm.dat, file="biclm.Rdata") } else { load("biclm.Rdata") } coef.biclm.dat <- summary(biclm.dat)$coef[-1, ] coef.biclm.dat <- coef.biclm.dat[order(coef.biclm.dat[, 1]), ] my.barplot.lm( rownames(coef.biclm.dat), coef.biclm.dat[, 1], err=coef.biclm.dat[, 2] * qnorm(0.975)) ## ---- ## glmnet X <- as.matrix(dat.lm)[, -1] y <- dat.lm[, 1] glmnet.dat <- glmnet(X, y) cri <- mapply(function(lambda, df) { beta <- as.numeric(coef(glmnet.dat, lambda)) resid <- drop(y - cbind(rep(1, nrow(X)), X) %*% beta) log(var(resid)) + log(nrow(X)) * df / nrow(X) }, glmnet.dat$lambda, glmnet.dat$df) best.i <- which.min(cri) coef.glmnet.dat <- drop(coef(glmnet.dat, glmnet.dat$lambda[best.i]))[-1] coef.glmnet.dat <- coef.glmnet.dat[coef.glmnet.dat != 0] coef.glmnet.dat <- sort(coef.glmnet.dat) my.barplot.lm(names(coef.glmnet.dat), coef.glmnet.dat) ## ---- ## msgps msgps.dat <- msgps(as.matrix(dat.lm)[, -1], dat.lm[, 1], penalty="alasso") coef.msgps.dat <- drop(coef(msgps.dat, tuning=msgps.dat$dfbic_result$tuning))[-1] coef.msgps.dat <- coef.msgps.dat[coef.msgps.dat != 0] coef.msgps.dat <- sort(coef.msgps.dat) my.barplot.lm(names(coef.msgps.dat), coef.msgps.dat) [/code]
It’s much easier to unsnedtard when you put it that way!