【R】重回帰分析による予測モデルを使って住宅の購入価格を予想する

備忘録

最近、研究の合間に確率・統計を勉強し直してる一環でKaggleに登録した。(学部生の時にもっとこういうので遊んでおけばよかったと後悔している)

Kaggleといえば機械学習入門用のタイタニックの生存者予測で有名で実際に私もRを使ってトライしている。Rのパッケージというのは本当に便利なものでちょっとデータをいじってあとはニューラルネットワークのパッケージ(neuralnet)を使っただけでそこそこ良いスコアが取れた(上位35パーセント)
こちらは生きるか死ぬかのClassification判別分析で今回の趣旨とは違うのでまた改めて別の記事で書きたいと思う。

今回は判別分析ではなくRegression回帰分析で構築した予測モデルを使って同じくKaggleの入門者用データセットにある住宅の購入価格の予測を行いたいと思う。回帰分析というのは予測値が今回の”住宅の価格”のように連続した値となるような予測モデルを構築すること。

例えば敷地面積が広ければ広いほど住宅の購入価格が高くなることはイメージし易いが既存のそれらのデータセットから回帰分析を行うと敷地面積(説明変数、X軸)に対する住宅の購入価格(目的変数、Y軸)が予想できる。
Rには1つまたは複数の説明変数(独立変数)から目的変数(従属変数)を求めるための単回帰・重回帰分析に用いられるパッケージがデフォルトで入っているので今回はそれを使う。

なお、勉強にあたって日本から確率統計の参考書を持ってきたのだがけっこう役に立ったので記載しておく。

確率・統計(理工系の数学入門コース7)

原因を探る統計学ー共分散構造分析入門

Kaggleの各Notebook、特にこれ

以下の手順に従って予測モデルを構築していく。

データの取り込みと前処理

まずはデータの取り込み。Kaggleの入門者用コンペティションにはだいたい予測モデル構築のためのTrain.csvとそのモデルを使って予測するための目的変数(この場合住宅価格)が欠落したTest.csvが入っている。
Train <- read.csv(“Train.csv”)
str(Train)

でTrainの構造(中身)を見てみると目的変数である住宅価格(SalePrice)を除いた変数の数はなんと80個。しかも英語が不動産用語っぽくて分かりづらい。今回は簡潔に済ませるため量で表すことができない質的変数や離散値などを除いた量的変数、そこからさらに住宅価格と相関の高いものをピックアップして重回帰分析による予測モデルを構築する。住宅価格は敷地面積が広ければ広いほど高いであろうという直感に従って以下の4つの量的変数をピックアップした。

TotalBsmtSF(地下フロア面積), GarageArea(車庫面積), GrLivArea(居住空間面積), 1stFlrSF(1階フロア面積)

HousePriceとの相関係数はそれぞれ0.61, 0.62, 0.70, 0.60となる。
SalePrice, TotalBsmtSF, GrarageArea, GrLivArea, X1stFlrSFで新たにデータセットを定義し、SalePriceを目的変数として回帰予測モデルを構築していきたい。

欠損値の処理

既存のデータの中には値が不明(NaN)のものがありどう欠損値を埋めるかは予測モデルの精度に影響を与えるが今回ピックアップした変数の中には欠損値はわずかな数しか含まれていないのでとりあえずその変数の平均値を入れておく。

外れ値の除去

欠損値無しのデータセットが整ったところでSalePriceとそれぞれの変数における散布図を見てみる。最小二乗法による回帰分析ではしばしば以下の散布図にも見られるような飛んでいる値(外れ値)が分析結果に影響を与える。なので面積のわりに不自然に安くなっているところ(赤線の右側)は除去することにする。(事故物件なのだろうか。。。)

回帰分析

外れ値を除去したところでデータセットの中身をもう少し具体的に見ていく。

回帰分析を行うにあたっては以下の4つの性質をデータセットの変数が保持しているかどうかが重要となる。

・Linearity(線形性)

・Normality(正規分散性)

・Homoscedasticity(同一分散性)

・Multicollinearity(多重共線性)

データセットはすでに説明変数と相関があるものをピックアップしているため線形性はクリアしていることにする。後者の2つはまた改めて取り扱うことにして今回はNormality(正規分散性)に着目する。まずは各変数の分布を見てみる。説明変数である住宅価格SalePriceの分布図をhist(“data”)lines(density(“data”))でプロットしたものを以下に示す。

アメリカの住宅購入事情は詳しくないが低価格帯に分布が集中している一方でいわゆる超高級住宅の価格は青天井であることがこのグラフから読み取ることができて、これらは私たちの直感と一致していることだと思う。
左右対称の形を取る正規分布よりも少し左によっていることが分かる。このような場合、最頻値(モード)は平均値や中央値よりも低くなり歪度(Skewness)が正であるという、らしい。歪度が見られる分布は一般的に対数を取ると正規分布上の形に近づくことが知られているので対数を取る。

正規分布の形に近づいた。回帰分析においてはこれを目的変数とする。
4つの説明変数についてもそれぞれの分布を見ていく。

ここで正規分布上になっているかどうかを確認するために、理論的な分布との適合度を調べるQ-Qプロットを使ってみる。Rではqqnorm(data), qqplot(data)で図示できる。他にもRでは適合度検定としてコルモゴロフ–スミルノフ検定がks.test()でデフォルトで使える。
4つの変数に対するQ-Qプロットの結果を以下に示す。

直線上に乗っているほど正規分布の形であるといえる。さきほどの歪度が見られる分布は下2つのように少し反った形となる。上2つは散布図でも分かるように値が0である場合が多々見られる。地下室と車庫がないケースである。この場合真数が0となってしまうため対数は取れない。0を除いた部分ではおおむね正規分布上に乗っていると判断して地下室面積と車庫面積においてはこのままの値を使用する。下2つに対して対数を取ったあとの分布とQ-Qプロットを以下に示す。どちらも正規分布の様相になったのでこれを回帰分析に用いる。

Rでは回帰分析にはlm()関数を用いる。分布形状を指定する場合はglm()関数でもいい。今回はガウシアン(正規)分布なのでどっちでもいいと思う。

回帰分析による予測モデルができたらあとはpredict()関数を使えばモデルに沿った目的変数を出力することができる。

予測モデルの構築と精度の確認

Testデータに対して構築した予測モデルを適用し目的変数である住宅価格が出力できればKaggleに提出できる。ちなみに本予測モデルで提出した時のRMSEスコアは0.22(0に近いほど正確)で4600位中4000位とかなり下の方からのスタートであった。トップ層は0.1近辺までの精度を出してきているのでまぁ勉強を継続しつつ気長にチャレンジしていきたい。
提出する前に自分の作った予測モデルの正確さを知りたければTrainデータを2:1くらいに分割してそれぞれTrainデータ、Testデータとして扱うと予測モデルがどれくらいの精度を持つのか確認できる。
今回は歪度を持つ分布に対する対数変換がどの程度予測モデルに影響を与えるのかTrainデータの2/3を使った予測モデルを使用して確認してみた。Trainデータの1/3をTestデータとした時の予測データ(対数変換あり・なし)と本来の値の分布の比較を以下に示す。

元々の分布の歪度が正で平均値がモードよりも高くなっているから対数変換無しのデータでは実際の値よりも少し高めに予測されているということ、だと思う。よく見てみると予測データ(赤)よりも実際のデータ(黒)の方が左右に裾が広がっていてとくに住宅価格が高い領域を十分に予測し切れていないことが分かる。今回は独立変数として各敷地面積のみを採用したけど本来のデータセットには80個もの変数が用意されているのでカテゴリカル変数や離散値も含めて住居価格により大きな影響を与える変数を使って回帰分析すると改善が見込まれるのだろう。

将来的にはニューラルネットワークとか組み合わせたSpace Mappingまでやれれば研究にも少し役に立つかもしれない、くらいのモチベーションで細々と続けていきたいと思います。

以下、今回のコード。train.csv, test.csv さえ作業ディレクトリにあればたぶんグラフの描画までやってくれると思う。もし今回の手順やコードに改善点等あれば教えてください。

今回のコード
Train <- read.csv("train.csv")
Test <- read.csv("test.csv")

#To pick up the numeric features
Train_M <- Train[,c("Id","SalePrice","TotalBsmtSF","GrLivArea","OverallQual","X1stFlrSF",
                    "GarageArea")]
Test_M <- Test[,c("Id","TotalBsmtSF","GrLivArea","OverallQual","X1stFlrSF",
                  "GarageArea")]

#NaN is replaced by the mean of each feature
Test_M$GarageArea[is.na(Test_M$GarageArea)] <-
  mean(Test_M$GarageArea[!is.na(Test_M$GarageArea)])
Test_M$TotalBsmtSF[is.na(Test_M$TotalBsmtSF)] <-
  mean(Test_M$TotalBsmtSF[!is.na(Test_M$TotalBsmtSF)])


#to show a graph of the sale price dependence on some features
png("SalePriceDependancy.png")
par(oma=c(0,0,4,0))
par(mfrow=c(2,2))
plot(Train$TotalBsmtSF,Train$SalePrice)
abline(v=4000,col="red") #Outlier
plot(Train$GarageArea,Train$SalePrice)
abline(v=1200,col="red") #Outlier
plot(Train$GrLivArea,Train$SalePrice)
abline(v=4500,col="red") #Outlier
plot(Train$X1stFlrSF,Train$SalePrice)
abline(v=2600,col="red") #Outlier
mtext(side = 3, line = 1, outer=T, text="SalePrice Vs. Four Features", cex=2)
dev.off()


##Outlier value is eliminated
Train_M <- Train_M[(Train_M$TotalBsmtSF <= 4000) & (Train_M$GarageArea<=1200) &
                     (Train_M$GrLivArea <= 4500) & (Train_M$X1stFlrSF <= 2600),]

#To show the distribution of some features if it´s normally distributed or not.
png("Distribution.png")
par(oma=c(0,0,4,0))
par(mfrow=c(2,2))
plot(density(Train_M$TotalBsmtSF))
plot(density(Train_M$GarageArea))
plot(density(Train_M$GrLivArea))
plot(density(Train_M$X1stFlrSF))
dev.off()

png("QQplot.png")
par(oma=c(0,0,4,0))
par(mfrow=c(2,2))
qqnorm(Train_M$TotalBsmtSF)
qqline(Train_M$TotalBsmtSF)

qqnorm(Train_M$GarageArea)
qqline(Train_M$GarageArea)

qqnorm(Train_M$GrLivArea)
qqline(Train_M$GrLivArea)

qqnorm(Train_M$X1stFlrSF)
qqline(Train_M$X1stFlrSF)

dev.off()

#GarageArea and TotalBsmtSF has zeros, so taking the log is impossible
SalePrice.log <- log10(Train_M$SalePrice)
GrLivArea.log <- log10(Train_M$GrLivArea)
X1stFlrSF.log <- log10(Train_M$X1stFlrSF)
#TotalBsmtSF.log <- log10(Train_M$TotalBsmtSF) #inappropriate
#GarageArea.log <- log10(Train_M$GarageArea) #inappropriate

#Logarithm features distribution showing it has normal distribution, where lm regression is available
png("LogQQplot.png")
par(oma=c(0,0,4,0))
par(mfrow=c(2,2))

plot(density(Train_M$GrLivArea.log))
plot(density(Train_M$X1stFlrSF.log))

qqnorm(Train_M$GrLivArea.log)
qqline(Train_M$GrLivArea.log)

qqnorm(Train_M$X1stFlrSF.log)
qqline(Train_M$X1stFlrSF.log)

dev.off()


#to make the new data frame of logged data and combine it with original train data
LogData1 <- data.frame("SalePrice.log" = SalePrice.log,"GrLivArea.log" = GrLivArea.log,
                       "X1stFlrSF.log" = X1stFlrSF.log )
Train_M <- cbind(Train_M,LogData1)


#To apply same preparation to test data
GrLivArea.log.Test <- log10(Test_M$GrLivArea)
X1stFlrSF.log.Test <- log10(Test_M$X1stFlrSF)
LogData2 <- data.frame("GrLivArea.log" = GrLivArea.log.Test,
                       "X1stFlrSF.log" = X1stFlrSF.log.Test )
Test_M <- cbind(Test_M,LogData2)


##Employing glm function for regression model
res <- glm(SalePrice.log~GrLivArea.log + X1stFlrSF.log +
              TotalBsmtSF + GarageArea, Train_M2,family = "gaussian")

#to show the summary result of this regression model
summary(res)


###Prediction#####
SalePrice.log.Test <- predict(res1,Test_M)
SalePrice.Test <- 10**SalePrice.log.Test


#####MakingNewData######
SalePrice_Test2 <- data.frame("SalePrice" = SalePrice.Test2, "SalePrice.log" = SalePrice.log.Test2)
SalePrice_Test3 <- data.frame("SalePrice" = SalePrice.Test3, "SalePrice.log" = SalePrice.log.Test3)
Test_M2 <- cbind(Test_M2,SalePrice_Test2)
Test_M3 <- cbind(Test_M3,SalePrice_Test3)
Test_New <- rbind(Test_M2,Test_M3)

#####MakingSubmission File#######
Submission <- data.frame("Id" = Test_M$Id, "SalePrice" = SalePrice.Test)
write.csv(Submission,"Submission_lm_New.csv",row.names=F)

コメント

タイトルとURLをコピーしました