Edited at

kaggleのkernels学習ノート~House prices: Lasso, XGBoost, and a detailed EDA~

題名の通りkaggleのkernelsを和訳したものです。自分の勉強の備忘録をまとめていきます。そして、データ分析に携わる誰かのお役に立てれば幸いです。和訳とか得意じゃないので間違っていたらごめんなさい。

今回はErik BruinさんのHouse prices: Lasso, XGBoost, and a detailed EDA

を和訳していきます。


1 Executive Summary

データセットをよく理解することに焦点を当て、このコンペティションに参加しました。EDAでは、多くの視覚化が含まれています。このバージョンにはモデリングも含まれています。


  • Lasso回帰は、RMSEスコアが0.1121と最もよくワークしました。変数の間にマルチコがあるということは、予想できました。Lasso回帰では、想定されているように、そのモデルで利用可能な変数をかなり限定します。

  • XGBoostモデルは、RMSEが0.1162と非常によくワークしました。

  • これら2つのアルゴリズムは非常に異なりますし、予測結果を平均化することで、予測を改善する可能性があります。

Kaggleはこのコンペについて以下のように説明しています。

自宅のバイヤーに自分の夢の家について説明するように依頼してください。おそらく地下の天井の高さや東西鉄道の近くから始まることはありません。 しかし、この競技場の競技会のデータセットは、寝室の数やホワイトピケットの柵の価格交渉よりもはるかに影響を与えることを証明しています。

アイオワ州のエイムズにある(ほぼ)住宅のあらゆる側面を説明する79の説明変数により、この競争は各家の最終価格を予測することに挑戦します。バリデーションのRMSEが、XGBoostのCVスコアより優れているので、私はLasso回帰の結果に2倍の重みを付けることにしました。


2 Introduction

Kaggleでは、このコンペについて以下のように説明されています。

不動産バイヤーに自分のドリームハウスについて説明するように尋ねてみてください。おそらく地下の天井の高さや東西にある鉄道の近くの話から始まることはありません。 しかし、このコンペのデータセットは、寝室の数やホワイトピケットの柵よりも、はるかに価格交渉に影響を与えることを証明しています。アイオワ州のエイムズにある(ほぼ)住宅のあらゆる側面を説明する79の説明変数により、このコンペでは、各家の最終価格を予測することに挑戦します。

スクリーンショット 0030-09-14 20.57.56.png


3 Loading and Exploring Data


3.1 Loading libraries required and reading the data into R

base以外に使用するRパッケージをロードします。


R

library(knitr)

library(ggplot2)
library(plyr)
library(dplyr)
library(corrplot)
library(caret)
library(gridExtra)
library(scales)
library(Rmisc)
library(ggrepel)
library(randomForest)
library(psych)
library(xgboost)

CSVファイルをRに読み込みます。


R

train <- read.csv("../input/train.csv", stringsAsFactors = F)

test <- read.csv("../input/test.csv", stringsAsFactors = F)


3.2 Data size and structure

訓練データは、文字型および整数型で構成されます。文字変数の大部分は実際には(序数的な)要素ですが、文字列としてRに読み込むことを選択しました。その大部分は、まずクリーニングやフィーチャエンジニアリングが必要です。合計で81個の列/変数があり、そのうちの最後の変数は応答変数(SalePrice)です。以下では、変数を垣間見ているだけです。それらのすべてについては、ドキュメントを通してより詳細に議論されています。


R

dim(train)

## [1] 1460 81
str(train[,c(1:10, 81)]) #display first 10 variables and the response variable
## 'data.frame': 1460 obs. of 11 variables:
## $ Id : int 1 2 3 4 5 6 7 8 9 10 ...
## $ MSSubClass : int 60 20 60 70 60 50 20 60 50 190 ...
## $ MSZoning : chr "RL" "RL" "RL" "RL" ...
## $ LotFrontage: int 65 80 68 60 84 85 75 NA 51 50 ...
## $ LotArea : int 8450 9600 11250 9550 14260 14115 10084 10382 6120 7420 ...
## $ Street : chr "Pave" "Pave" "Pave" "Pave" ...
## $ Alley : chr NA NA NA NA ...
## $ LotShape : chr "Reg" "Reg" "IR1" "IR1" ...
## $ LandContour: chr "Lvl" "Lvl" "Lvl" "Lvl" ...
## $ Utilities : chr "AllPub" "AllPub" "AllPub" "AllPub" ...
## $ SalePrice : int 208500 181500 223500 140000 250000 143000 307000 200000 129900 118000 ...

#Getting rid of the IDs but keeping the test IDs in a vector. These are needed to compose the submission file

test_labels <- test$Id
test$Id <- NULL
train$Id <- NULL
test$SalePrice <- NA

all <- rbind(train, test)
dim(all)
## [1] 2919 80


IDがなければ、データフレームは79個の予測変数と応答変数(SalePrice)で構成されます。


4 Exploring some of the most important variables


4.1 The response variable; SalePrice

見ての通り、販売価格は右に偏っています。これは、非常に高価な住宅を買う人がほとんどいないと予想できます。これを念頭に置いて、モデル化する前に対処する必要があります。


R

ggplot(data=all[!is.na(all$SalePrice),], aes(x=SalePrice)) +

geom_histogram(fill="blue", binwidth = 10000) +
scale_x_continuous(breaks= seq(0, 800000, by=100000), labels = comma)

unnamed-chunk-6-1.png


R

summary(all$SalePrice)

## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 34900 129975 163000 180921 214000 755000 1459


4.2 The most important numeric predictors


4.2.1 Correlations with SalePrice

全体として、SalePriceとの相関が少なくとも0.5である数値変数が10個あります。これらの相関はすべて正の相関です。


R

numericVars <- which(sapply(all, is.numeric)) #index vector numeric variables

numericVarNames <- names(numericVars) #saving names vector for use later on
cat('There are', length(numericVars), 'numeric variables')
## There are 37 numeric variables
all_numVar <- all[, numericVars]
cor_numVar <- cor(all_numVar, use="pairwise.complete.obs")
#correlations of all numeric variables

#sort on decreasing correlations with SalePrice
cor_sorted <- as.matrix(sort(cor_numVar[,'SalePrice'], decreasing = TRUE))
#select only high corelations
CorHigh <- names(which(apply(cor_sorted, 1, function(x) abs(x)>0.5)))
cor_numVar <- cor_numVar[CorHigh, CorHigh]

corrplot.mixed(cor_numVar, tl.col="black", tl.pos = "lt")


unnamed-chunk-8-1.png

このセクションの残りの部分では、SalePriceとSalePriceとの相関が最も高い2つ(Overall QualityとAbove Grade Ground Living Area square feet)の予測子の間の関係を視覚化します。

マルチコが問題であることも明らかです。たとえば、GarageCarsとGarageAreaの間の相関は非常に高く(0.89)、どちらもSalePriceと同様に(高い)相関を持っています。SalePriceとの相関が0.5より高い6つの6つの変数は次のとおりです。


  • TotalBsmtSF:地下1階の合計平方フィート

  • 1stFlrSF:1階平方フィート

  • フルバス:グレード上のフルバスルーム

  • TotRmsAbvGrd:グレードを超える客室の合計(バスルームは含まれません )

  • YearBuilt:建設年

  • YearRemodAdd:リフォーム日


4.2.2 Overall Quality

SalesPriceとの相関が2番目に高い数値変数は、Above Grade Living Areaです。 これは多くの意味があります。大きな不動産は一般的に高価です。


R

ggplot(data=all[!is.na(all$SalePrice),], aes(x=GrLivArea, y=SalePrice))+

geom_point(col='blue') + geom_smooth(method = "lm", se=FALSE, color="black", aes(group=1)) +
scale_y_continuous(breaks= seq(0, 800000, by=100000), labels = comma) +
geom_text_repel(aes(label = ifelse(all$GrLivArea[!is.na(all$SalePrice)]>4500, rownames(all), '')))

unnamed-chunk-10-1.png

低いSalePricesを持つ2つの家は外れ値に見えます(524と1299の不動産。グラフのラベルを参照)。外れ値を取ることは危険であるかもしれないので、まだ取り除かないでおきます。たとえば、Overall Qualityに関する低いスコアは、低価格を説明することができます。しかし、以下に示すように、これらの2つの家は、実際にはOverall Qualityで最高点を獲得します。したがって、1299と524を外れ値として取り出す主要な候補として念頭に置いておきます。


4.2.3 Above Grade (Ground) Living Area (square feet)

SalesPriceとの相関が2番目に高い数値変数は、Above Grade Living Areaです。大きな家は一般的に高価です。


R

ggplot(data=all[!is.na(all$SalePrice),], aes(x=GrLivArea, y=SalePrice))+

geom_point(col='blue') + geom_smooth(method = "lm", se=FALSE, color="black", aes(group=1)) +
scale_y_continuous(breaks= seq(0, 800000, by=100000), labels = comma) +
geom_text_repel(aes(label = ifelse(all$GrLivArea[!is.na(all$SalePrice)]>4500, rownames(all), '')))

unnamed-chunk-10-1.png


R

all[c(524, 1299), c('SalePrice', 'GrLivArea', 'OverallQual')]

## SalePrice GrLivArea OverallQual
## 524 184750 4676 10
## 1299 160000 5642 10


5 Missing data, label encoding, and factorizing variables


5.1 Completeness of the data

まず、どの変数に欠損値が含まれているかを確認したいと思います。


R

NAcol <- which(colSums(is.na(all)) > 0)

sort(colSums(sapply(all[NAcol], is.na)), decreasing = TRUE)
## PoolQC MiscFeature Alley Fence SalePrice
## 2909 2814 2721 2348 1459
## FireplaceQu LotFrontage GarageYrBlt GarageFinish GarageQual
## 1420 486 159 159 159
## GarageCond GarageType BsmtCond BsmtExposure BsmtQual
## 159 157 82 82 81
## BsmtFinType2 BsmtFinType1 MasVnrType MasVnrArea MSZoning
## 80 79 24 23 4
## Utilities BsmtFullBath BsmtHalfBath Functional Exterior1st
## 2 2 2 2 1
## Exterior2nd BsmtFinSF1 BsmtFinSF2 BsmtUnfSF TotalBsmtSF
## 1 1 1 1 1
## Electrical KitchenQual GarageCars GarageArea SaleType
## 1 1 1 1 1

cat('There are', length(NAcol), 'columns with missing values')
## There are 35 columns with missing values


SalePriceの1459のNAsはテストデータのサイズと完全に一致します。また、34個の予測変数でNAを修正しなければなりません。


5.2 Imputing missing data

このセクションでは、欠損値を含む34個の予測変数を修正します。他の変数とグループを形成できる変数を見つけたら、それらをグループとして扱います。たとえば、プール、ガレージ、および地下に関連する複数の変数があります。

文書をできるだけ読みやすくするために、knitrが提供する「タブ」オプションを使用することにしました。それぞれのタブの下にある変数について、簡単な分析を見ることができます。すべてのセクションを見る必要はなく、いくつかのタブを見ることもできます。もしそうなら、私は特に、garagesとBasementのセクションは価値があると思います。

文字型も明確な順序があれば順序付き変数に、順序がないカテゴリであればファクタ型に変換します。後でone-hotエンコーディング(model.matrix関数を使用)を使って、これらのファクタ型を数値に変換します。

※ここではgaragesのセクションを和訳しておきます。それ以外にセクションに興味がある方は、リンク先から確認してください。


5.2.7 Altogether, there are 7 variables related to garages

そのうちの2つにはNA(GarageCarsとGarageArea)が1つ、157のNAs(GarageType)があり、4つの変数には159のNAがあります。まず最初にGarageYrBltの159件の欠損すべて置き換える予定です。Year garageはYearBuiltの値で構築されました。


R

all$GarageYrBlt[is.na(all$GarageYrBlt)] <- all$YearBuilt[is.na(all$GarageYrBlt)]


NAsは「ガレージなし」を意味するので、他の3つの変数との差異がどこから来るのかを知りたい。


R

#check if all 157 NAs are the same observations among the variables with 157/159 NAs

length(which(is.na(all$GarageType) & is.na(all$GarageFinish) & is.na(all$GarageCond) & is.na(all$GarageQual)))
## [1] 157
#Find the 2 additional NAs
kable(all[!is.na(all$GarageType) & is.na(all$GarageFinish), c('GarageCars', 'GarageArea', 'GarageType', 'GarageCond', 'GarageQual', 'GarageFinish')])

GarageCars
GarageArea
GarageType
GarageCond
GarageQual
GarageFinish

2127
1
360
Detchd
NA
NA
NA

2577
NA
NA
Detchd
NA
NA
NA

GarageType内の157のNAsは、GarageCondition、GarageQuality、およびGarageFinishでもNAになります。 違いは、不動産2127と不動産2577にあります。ご覧のとおり、2127は実際にガレージを持っているようですが、不動産2577はありません。したがって、ガレージなしの158個の不動産があるはずです。不動産2127を修正するために、GarageCond、GarageQual、およびGarageFinishの最も一般的な値(モード)を代入します。


R

#Imputing modes.

all$GarageCond[2127] <- names(sort(-table(all$GarageCond)))[1]
all$GarageQual[2127] <- names(sort(-table(all$GarageQual)))[1]
all$GarageFinish[2127] <- names(sort(-table(all$GarageFinish)))[1]

#display "fixed" house
kable(all[2127, c('GarageYrBlt', 'GarageCars', 'GarageArea', 'GarageType', 'GarageCond', 'GarageQual', 'GarageFinish')])


GarageYrBlt
GarageCars
GarageArea
GarageType
GarageCond
GarageQual
GarageFinish

2127
1910
1
360
Detchd
TA
TA
Unf


GarageCarsとGarageArea

どちらも1つのNAです。上記のように、不動産2577がこれにあたります。GarageTypeが "detached"で、他のすべてのGarage変数がこの家屋にGarageがないことを示しているように見えるため、問題が発生した可能性があります。


R

#fixing 3 values for house 2577

all$GarageCars[2577] <- 0
all$GarageArea[2577] <- 0
all$GarageType[2577] <- NA

#check if NAs of the character variables are now all 158
length(which(is.na(all$GarageType) & is.na(all$GarageFinish) & is.na(all$GarageCond) & is.na(all$GarageQual)))
## [1] 158


今、ガレージに関連する4つの変数はすべて、「No Garage」に対応する同じ158個のNAsを持っています。このセクションの残りの部分ですべてを修正します。


GarageType: Garage location

値は序数ではないように見えるので、ファクタに変換します。


R

2Types   More than one type of garage

Attchd Attached to home
Basment Basement Garage
BuiltIn Built-In (Garage part of house - typically has room above garage)
CarPort Car Port
Detchd Detached from home
NA No Garage
all$GarageType[is.na(all$GarageType)] <- 'No Garage'
all$GarageType <- as.factor(all$GarageType)
table(all$GarageType)
##
## 2Types Attchd Basment BuiltIn CarPort Detchd No Garage
## 23 1723 36 186 15 778 158


GarageFinish: Interior finish of the garage

順序変数にします。


R

 Fin  Finished

RFn Rough Finished
Unf Unfinished
NA No Garage

all$GarageFinish[is.na(all$GarageFinish)] <- 'None'
Finish <- c('None'=0, 'Unf'=1, 'RFn'=2, 'Fin'=3)

all$GarageFinish<-as.integer(revalue(all$GarageFinish, Finish))
table(all$GarageFinish)
##
## 0 1 2 3
## 158 1231 811 719



GarageQual: Garage quality

品質のベクトルとして、順序付けしておきます。


R

Ex   Excellent

Gd Good
TA Typical/Average
Fa Fair
Po Poor
NA No Garage

all$GarageQual[is.na(all$GarageQual)] <- 'None'
all$GarageQual<-as.integer(revalue(all$GarageQual, Qualities))
table(all$GarageQual)
##
## 0 1 2 3 4 5
## 158 5 124 2605 24 3



GarageCond: Garage condition

品質のベクトルとして、順序付けしておきます。


R

Ex   Excellent

Gd Good
TA Typical/Average
Fa Fair
Po Poor
NA No Garage

all$GarageCond[is.na(all$GarageCond)] <- 'None'
all$GarageCond<-as.integer(revalue(all$GarageCond, Qualities))
table(all$GarageCond)
##
## 0 1 2 3 4 5
## 158 14 74 2655 15 3



5.3 Label encoding/factorizing the remaining character variables

この時点で、NAを持つすべての変数が処理されていることを確認しました。しかし、まだ欠損値のない残りの文字型変数を処理する必要があります。前のセクションと同様に、私は変数のグループのためのタブを作成しました。


R

Charcol <- names(all[,sapply(all, is.character)])

Charcol
## [1] "Street" "LandContour" "LandSlope" "Neighborhood"
## [5] "Condition1" "Condition2" "BldgType" "HouseStyle"
## [9] "RoofStyle" "RoofMatl" "Foundation" "Heating"
## [13] "HeatingQC" "CentralAir" "PavedDrive"

cat('There are', length(Charcol), 'remaining columns with character values')
## There are 15 remaining columns with character values


ここでは、foundationのセクションを和訳しておきます。


5.3.1 Foundation: foundationの種類


R

BrkTil          Brick & Tile

CBlock Cinder Block
PConc Poured Contrete
Slab Slab
Stone Stone
Wood Wood

#No ordinality, so converting into factors
all$Foundation <- as.factor(all$Foundation)
table(all$Foundation)
##
## BrkTil CBlock PConc Slab Stone Wood
## 311 1235 1308 49 11 5
sum(table(all$Foundation))
## [1] 2919



5.4 Changing some numeric variables into factors

この時点で、すべての変数は完全(No NAs)であり、すべての文字型変数は因子の数値ラベルに変換されます。しかし、数値を記録する3つの変数がありますが、実際にはカテゴリに分類する必要があります。


5.4.1 Year and Month Sold

YearBuilt内の順序化は理にかなっていますが(古い家はそれほど価値がありません)、ここでは5年間の販売について話しています。これらの年には経済危機も含まれています。たとえば、2009年の販売価格(崩壊後)は、2007年よりもはるかに低い可能性が非常に高いです。YrSoldをモデリングする前に、ファクタ型に変換する必要がありますが、Age変数を作成するにはYrSoldの数値版が必要です。まだそれはやっていません。Month Soldも整数の変数です。 しかし、12月は1月より「良い」わけではありません。したがって、私はMoSold値をファクタ型に戻します。


R

str(all$YrSold)

## int [1:2919] 2008 2007 2008 2006 2008 2009 2007 2009 2008 2008 ...

str(all$MoSold)
## int [1:2919] 2 5 9 2 12 10 8 11 4 1 ...
all$MoSold <- as.factor(all$MoSold)


予想より少し緩やかですが、2007年末の銀行危機の影響は確かに見られます。2007年を境に、価格は徐々に低下しています。しかし、以下のように、季節性はより大きな役割を果たしているようです。


R

ys <- ggplot(all[!is.na(all$SalePrice),], aes(x=as.factor(YrSold), y=SalePrice)) +

geom_bar(stat='summary', fun.y = "median", fill='blue')+
scale_y_continuous(breaks= seq(0, 800000, by=25000), labels = comma) +
geom_label(stat = "count", aes(label = ..count.., y = ..count..)) +
coord_cartesian(ylim = c(0, 200000)) +
geom_hline(yintercept=163000, linetype="dashed", color = "red") #dashed line is median SalePrice

ms <- ggplot(all[!is.na(all$SalePrice),], aes(x=MoSold, y=SalePrice)) +
geom_bar(stat='summary', fun.y = "median", fill='blue')+
scale_y_continuous(breaks= seq(0, 800000, by=25000), labels = comma) +
geom_label(stat = "count", aes(label = ..count.., y = ..count..)) +
coord_cartesian(ylim = c(0, 200000)) +
geom_hline(yintercept=163000, linetype="dashed", color = "red") #dashed line is median SalePrice

grid.arrange(ys, ms, widths=c(1,2))


unnamed-chunk-85-1.png


5.4.2 MSSubClass

MSSubClass:販売に関与する住居のタイプを識別します。


R

  20  1-STORY 1946 & NEWER ALL STYLES

30 1-STORY 1945 & OLDER
40 1-STORY W/FINISHED ATTIC ALL AGES
45 1-1/2 STORY - UNFINISHED ALL AGES
50 1-1/2 STORY FINISHED ALL AGES
60 2-STORY 1946 & NEWER
70 2-STORY 1945 & OLDER
75 2-1/2 STORY ALL AGES
80 SPLIT OR MULTI-LEVEL
85 SPLIT FOYER
90 DUPLEX - ALL STYLES AND AGES
120 1-STORY PUD (Planned Unit Development) - 1946 & NEWER
150 1-1/2 STORY PUD - ALL AGES
160 2-STORY PUD - 1946 & NEWER
180 PUD - MULTILEVEL - INCL SPLIT LEV/FOYER
190 2 FAMILY CONVERSION - ALL STYLES AND AGES

これらのクラスは数値としてコード化されていますが、実際はカテゴリです。


R

str(all$MSSubClass)

## int [1:2919] 60 20 60 70 60 50 20 60 50 190 ...
all$MSSubClass <- as.factor(all$MSSubClass)

#revalue for better readability
all$MSSubClass<-revalue(all$MSSubClass, c('20'='1 story 1946+', '30'='1 story 1945-', '40'='1 story unf attic', '45'='1,5 story unf', '50'='1,5 story fin', '60'='2 story 1946+', '70'='2 story 1945-', '75'='2,5 story all ages', '80'='split/multi level', '85'='split foyer', '90'='duplex all style/age', '120'='1 story PUD 1946+', '150'='1,5 story PUD all', '160'='2 story PUD 1946+', '180'='PUD multilevel', '190'='2 family conversion'))

str(all$MSSubClass)
## Factor w/ 16 levels "1 story 1946+",..: 6 1 6 7 6 5 1 6 5 16 ...



6 Visualization of important variables

私は最終的に、すべての文字型変数がファクタ型に変換されたか、ラベルが数字にコード化するところまで処理しました。さらに、3つの数値変数がファクタに変換され、1つの変数が削除されました。以下に示すように、数値変数の数は56(応答変数を含む)になり、残りの23の変数はカテゴリになります。


R

numericVars <- which(sapply(all, is.numeric)) #index vector numeric variables

factorVars <- which(sapply(all, is.factor)) #index vector factor variables
cat('There are', length(numericVars), 'numeric variables, and', length(factorVars), 'categoric variables')

## There are 56 numeric variables, and 23 categoric variables



6.1 Correlations again

相関関係をもう一度チェックしています。 ご覧のとおり、SalePriceとの相関が0.5以上の変数の数は、10(セクション4.2.1を参照)から16に増加しました。


R

all_numVar <- all[, numericVars]

cor_numVar <- cor(all_numVar, use="pairwise.complete.obs") #correlations of all numeric variables

#sort on decreasing correlations with SalePrice
cor_sorted <- as.matrix(sort(cor_numVar[,'SalePrice'], decreasing = TRUE))
#select only high corelations
CorHigh <- names(which(apply(cor_sorted, 1, function(x) abs(x)>0.5)))
cor_numVar <- cor_numVar[CorHigh, CorHigh]

corrplot.mixed(cor_numVar, tl.col="black", tl.pos = "lt", tl.cex = 0.7,cl.cex = .7, number.cex=.7)


unnamed-chunk-88-1.png


6.2 Finding variable importance with a quick Random Forest

相関関係は、最も重要な数値変数とそれらの変数間の概要を教えてくれますが、視覚化に移行する前にカテゴリ変数を含む最も重要な変数の概要を知りたいと考えました。

私は、パッケージのcalc.relimp関数を使って素早く線形回帰モデルを使って変数の相対的重要性を得ようとしました。また、変数を重要であるかどうかを分けるパッケージborutaboruta関数を試しました。

しかし、これらの方法では多くの時間がかかりました。私は可変的な重要性の兆候を得たいだけなので、結局はそれを単純にして、素早くランダムフォレストモデルを使用することに決めました。


R

set.seed(2018)

quick_RF <- randomForest(x=all[1:1460,-79], y=all$SalePrice[1:1460], ntree=100,importance=TRUE)
imp_RF <- importance(quick_RF)
imp_DF <- data.frame(Variables = row.names(imp_RF), MSE = imp_RF[,1])
imp_DF <- imp_DF[order(imp_DF$MSE, decreasing = TRUE),]

ggplot(imp_DF[1:20,], aes(x=reorder(Variables, MSE), y=MSE, fill=MSE)) + geom_bar(stat = 'identity') + labs(x = 'Variables', y= '% increase MSE if variable is randomly permuted') + coord_flip() + theme(legend.position="none")


unnamed-chunk-89-1.png


6.2.1 Ground Living Areaとother surface related variables

最初の探索では、Ground Living AreaとSalePriceの関係はすでに視覚化しているので、今は分布自体を表示します。トップ20にはさらに平方フィートの表面測定値があるので、私はこのセクションでそれらをまとめる機会を得ています。GarageAreaはGarageの変数セクションで処理されている点は注意してください。この変数は地上居住地(0.81)と高い相関があるため、TotRmsAbvGrdも追加しています。


R

s1 <- ggplot(data= all, aes(x=GrLivArea)) +

geom_density() + labs(x='Square feet living area')
s2 <- ggplot(data=all, aes(x=as.factor(TotRmsAbvGrd))) +
geom_histogram(stat='count') + labs(x='Rooms above Ground')
s3 <- ggplot(data= all, aes(x=X1stFlrSF)) +
geom_density() + labs(x='Square feet first floor')
s4 <- ggplot(data= all, aes(x=X2ndFlrSF)) +
geom_density() + labs(x='Square feet second floor')
s5 <- ggplot(data= all, aes(x=TotalBsmtSF)) +
geom_density() + labs(x='Square feet basement')
s6 <- ggplot(data= all[all$LotArea<100000,], aes(x=LotArea)) +
geom_density() + labs(x='Square feet lot')
s7 <- ggplot(data= all, aes(x=LotFrontage)) +
geom_density() + labs(x='Linear feet lot frontage')
s8 <- ggplot(data= all, aes(x=LowQualFinSF)) +
geom_histogram() + labs(x='Low quality square feet 1st & 2nd')

layout <- matrix(c(1,2,5,3,4,8,6,7),4,2,byrow=TRUE)
multiplot(s1, s2, s3, s4, s5, s6, s7, s8, layout=layout)


unnamed-chunk-90-1.png

外れ値のためのこれらの変数のいくつかを調査します。ロットの視覚化のために、私はすでに100,000平方フィート(4軒)以上のロットを取り出しました。

GrLivAreaは1階と2階の平方フィートのちょうど合計のようでした。しかし、それ以降のバージョンでは、LowQualFinSFと呼ばれる変数もあることがわかりました。

上記のように(低品質の平方フィートの1階と2階)ほぼすべての家にはこれがありません(わずか40の家にはいくつかあります)。これらの平方フィートは実際にGrLivAreaに含まれていることが分かります。 これらの3つの変数とGrLivAreaの間の相関は厳密に1となります。


R

cor(all$GrLivArea, (all$X1stFlrSF + all$X2ndFlrSF + all$LowQualFinSF))

## [1] 1
head(all[all$LowQualFinSF>0, c('GrLivArea', 'X1stFlrSF', 'X2ndFlrSF', 'LowQualFinSF')])
## GrLivArea X1stFlrSF X2ndFlrSF LowQualFinSF
## 52 1176 816 0 360
## 89 1526 1013 0 513
## 126 754 520 0 234
## 171 1382 854 0 528
## 186 3608 1518 1518 572
## 188 1656 808 704 144


6.2.2 The most important categorical variable; Neighborhood

最初のグラフは、NeighorhoodによるSalePriceの中央値を示しています。 訓練データ内の各近所の頻度(家の数)がラベルに表示されます。下の2番目のグラフは、すべてのデータの頻度を示しています。


R

n1 <- ggplot(all[!is.na(all$SalePrice),], aes(x=Neighborhood, y=SalePrice)) +

geom_bar(stat='summary', fun.y = "median", fill='blue') +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_y_continuous(breaks= seq(0, 800000, by=50000), labels = comma) +
geom_label(stat = "count", aes(label = ..count.., y = ..count..), size=3) +
geom_hline(yintercept=163000, linetype="dashed", color = "red") #dashed line is median SalePrice
n2 <- ggplot(data=all, aes(x=Neighborhood)) +
geom_histogram(stat='count')+
geom_label(stat = "count", aes(label = ..count.., y = ..count..), size=3)+
theme(axis.text.x = element_text(angle = 45, hjust = 1))
grid.arrange(n1, n2)

unnamed-chunk-92-1.png


6.2.3 Overall Quality, and other Quality variables

初期の調査で既にOverall QualityとSalePriceの関係を視覚化していますが、頻度分布も視覚化したいと思います。より多くの品質測定値があるので、私はこのセクションでそれらを要約しておきます。


R

q1 <- ggplot(data=all, aes(x=as.factor(OverallQual))) +

geom_histogram(stat='count')
q2 <- ggplot(data=all, aes(x=as.factor(ExterQual))) +
geom_histogram(stat='count')
q3 <- ggplot(data=all, aes(x=as.factor(BsmtQual))) +
geom_histogram(stat='count')
q4 <- ggplot(data=all, aes(x=as.factor(KitchenQual))) +
geom_histogram(stat='count')
q5 <- ggplot(data=all, aes(x=as.factor(GarageQual))) +
geom_histogram(stat='count')
q6 <- ggplot(data=all, aes(x=as.factor(FireplaceQu))) +
geom_histogram(stat='count')
q7 <- ggplot(data=all, aes(x=as.factor(PoolQC))) +
geom_histogram(stat='count')

layout <- matrix(c(1,2,8,3,4,8,5,6,7),3,3,byrow=TRUE)
multiplot(q1, q2, q3, q4, q5, q6, q7, layout=layout)


unnamed-chunk-93-1.png

Overall Qualityは非常に重要であり、他の変数よりもきめ細かなものです。External Qualityも重要ですが、Overall Quality(0.73)と高い相関があります。Kitchen Qualityは、すべての家屋にキッチンがあり、いくつかのわずかに違います。Garage Qualityは、ガレージの大半がQ3を持っているので、それほど違いはないようです。Fireplace Qualityは、高い相関関係のリストにあり、重要な変数リストに含まれています。PoolQCは非常にスパースです。 後でプールを持つ変数を作成します。


6.2.4 The second most important categorical variable; MSSubClass

最初のプロットは、MSSubClassによるSalePriceの中央値を示します。訓練データ内の各MSSubClassの頻度(住居数)がラベルに表示されます。ヒストグラムはすべてのデータの頻度を示します。ほとんどの家は比較的新しいもので、1〜2階建てです。


R

ms1 <- ggplot(all[!is.na(all$SalePrice),], aes(x=MSSubClass, y=SalePrice)) +

geom_bar(stat='summary', fun.y = "median", fill='blue') +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_y_continuous(breaks= seq(0, 800000, by=50000), labels = comma) +
geom_label(stat = "count", aes(label = ..count.., y = ..count..), size=3) +
geom_hline(yintercept=163000, linetype="dashed", color = "red") #dashed line is median SalePrice
ms2 <- ggplot(data=all, aes(x=MSSubClass)) +
geom_histogram(stat='count')+
geom_label(stat = "count", aes(label = ..count.., y = ..count..), size=3) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
grid.arrange(ms1, ms2)

unnamed-chunk-94-1.png


6.2.5 Garage variables

Garageのいくつかの変数はSalePriceと高い相関性を持ち、ランダムフォレストの上位20のリストにもあります。しかし、それらの間にはマルチコがあり、7つのガレージ変数がとにかく多すぎると思います。3つの変数のようなもので十分です。しかし、私が選択する前に、私はこのセクションでそれらのすべてを視覚化しています。


R

#correct error

all$GarageYrBlt[2593] <- 2007 #this must have been a typo. GarageYrBlt=2207, YearBuilt=2006, YearRemodAdd=2007.
g1 <- ggplot(data=all[all$GarageCars !=0,], aes(x=GarageYrBlt)) +
geom_histogram()
g2 <- ggplot(data=all, aes(x=as.factor(GarageCars))) +
geom_histogram(stat='count')
g3 <- ggplot(data= all, aes(x=GarageArea)) +
geom_density()
g4 <- ggplot(data=all, aes(x=as.factor(GarageCond))) +
geom_histogram(stat='count')
g5 <- ggplot(data=all, aes(x=GarageType)) +
geom_histogram(stat='count')
g6 <- ggplot(data=all, aes(x=as.factor(GarageQual))) +
geom_histogram(stat='count')
g7 <- ggplot(data=all, aes(x=as.factor(GarageFinish))) +
geom_histogram(stat='count')

layout <- matrix(c(1,5,5,2,3,8,6,4,7),3,3,byrow=TRUE)
multiplot(g1, g2, g3, g4, g5, g6, g7, layout=layout)


unnamed-chunk-96-1.png

セクション4.2ですでに述べたように、GarageCarsとGarageAreaは高く相関しています。 ここでは、GarageQualとGarageCondも高い相関を示しており、両方ともレベル=3が多くを占めています。


6.2.6 Basement variables

Garage変数と同様に、複数のBasement変数が相関行列とTop20リストで重要なようです。しかし、11のBasement変数は過度のようです。何かを決める前に、それらの8つを下に視覚化しています。2つのBathroom変数は、フィーチャエンジニアリング(7.1節)で扱われ、 Basement square feetは6.2.1節ですでに説明されています。


R

b1 <- ggplot(data=all, aes(x=BsmtFinSF1)) +

geom_histogram() + labs(x='Type 1 finished square feet')
b2 <- ggplot(data=all, aes(x=BsmtFinSF2)) +
geom_histogram()+ labs(x='Type 2 finished square feet')
b3 <- ggplot(data=all, aes(x=BsmtUnfSF)) +
geom_histogram()+ labs(x='Unfinished square feet')
b4 <- ggplot(data=all, aes(x=as.factor(BsmtFinType1))) +
geom_histogram(stat='count')+ labs(x='Rating of Type 1 finished area')
b5 <- ggplot(data=all, aes(x=as.factor(BsmtFinType2))) +
geom_histogram(stat='count')+ labs(x='Rating of Type 2 finished area')
b6 <- ggplot(data=all, aes(x=as.factor(BsmtQual))) +
geom_histogram(stat='count')+ labs(x='Height of the basement')
b7 <- ggplot(data=all, aes(x=as.factor(BsmtCond))) +
geom_histogram(stat='count')+ labs(x='Rating of general condition')
b8 <- ggplot(data=all, aes(x=as.factor(BsmtExposure))) +
geom_histogram(stat='count')+ labs(x='Walkout or garden level walls')

layout <- matrix(c(1,2,3,4,5,9,6,7,8),3,3,byrow=TRUE)
multiplot(b1, b2, b3, b4, b5, b6, b7, b8, layout=layout)


unnamed-chunk-97-1.png


7 Feature engineering


7.1 Total number of Bathrooms

ここに4つのBathroom変数があります。個別では、これらの変数はあまり重要ではありません。しかし、それらを1つの予測子にすると、これは強力なものになる可能性が高いと仮定します。「パウダールームやゲストバスとしても知られているハーフバスは、4つの主要なバスルームコンポーネントのうち、トイレとシンクの2つしかありません」。したがって、私はまた、半分のバスルームを半分とします。


R

all$TotBathrooms <- all$FullBath + (all$HalfBath*0.5) + all$BsmtFullBath + (all$BsmtHalfBath*0.5)


最初のグラフで分かるように、明らかに相関があるようです(0.63)。2番目のグラフには、すべてのデータのBathroomsの頻度分布が示されています。


R

tb1 <- ggplot(data=all[!is.na(all$SalePrice),], aes(x=as.factor(TotBathrooms), y=SalePrice))+

geom_point(col='blue') + geom_smooth(method = "lm", se=FALSE, color="black", aes(group=1)) +
scale_y_continuous(breaks= seq(0, 800000, by=100000), labels = comma)
tb2 <- ggplot(data=all, aes(x=as.factor(TotBathrooms))) +
geom_histogram(stat='count')
grid.arrange(tb1, tb2)

unnamed-chunk-99-1.png


7.2 Adding ‘House Age’, ‘Remodeled (Yes/No)’, and IsNew variables

全体として、House Ageに関係する3つの変数があります。YearBlt、YearRemodAdd、およびYearSoldです。Remodeling/Additionがない場合、YearRemodAddのデフォルトはYearBuiltになります。AgeRemodeledとYearSoldを使用して年齢を決定します。しかし、古い建築物の一部は常に残っており、家の一部のみが改装されている可能性があるため、改装済みのはい/いいえ変数も紹介します。 これは、年齢がリモデリングの日付に基づいている場合、同じ年の最初から作成された家よりも価値が低いことを示す、何らかの種類のペナルティパラメータと見なされます。


R

all$Remod <- ifelse(all$YearBuilt==all$YearRemodAdd, 0, 1) #0=No Remodeling, 1=Remodeling

all$Age <- as.numeric(all$YrSold)-all$YearRemodAdd
ggplot(data=all[!is.na(all$SalePrice),], aes(x=Age, y=SalePrice))+
geom_point(col='blue') + geom_smooth(method = "lm", se=FALSE, color="black", aes(group=1)) +
scale_y_continuous(breaks= seq(0, 800000, by=100000), labels = comma)

unnamed-chunk-101-1.png

予想どおり、グラフは年齢との負の相関を示しています(古い家はそれほど価値がありません)。


R

cor(all$SalePrice[!is.na(all$SalePrice)], all$Age[!is.na(all$SalePrice)])

## [1] -0.5090787

以下に見られるように、期待通り、改装された家屋は価値がありません。


R

ggplot(all[!is.na(all$SalePrice),], aes(x=as.factor(Remod), y=SalePrice)) +

geom_bar(stat='summary', fun.y = "median", fill='blue') +
geom_label(stat = "count", aes(label = ..count.., y = ..count..), size=6) +
scale_y_continuous(breaks= seq(0, 800000, by=50000), labels = comma) +
theme_grey(base_size = 18) +
geom_hline(yintercept=163000, linetype="dashed") #dashed line is median SalePrice

unnamed-chunk-103-1.png

最後に、IsNew変数を作成しています。 全体として、データセットには116の新しい家があります。


R

all$IsNew <- ifelse(all$YrSold==all$YearBuilt, 1, 0)

table(all$IsNew)
##
## 0 1
## 2803 116

これらの116の新しい家は訓練データとテストデータの間にかなり均等に配分されています。新しい家はかなり平均的に価値があります。


R

ggplot(all[!is.na(all$SalePrice),], aes(x=as.factor(IsNew), y=SalePrice)) +

geom_bar(stat='summary', fun.y = "median", fill='blue') +
geom_label(stat = "count", aes(label = ..count.., y = ..count..), size=6) +
scale_y_continuous(breaks= seq(0, 800000, by=50000), labels = comma) +
theme_grey(base_size = 18) +
geom_hline(yintercept=163000, linetype="dashed") #dashed line is median SalePrice

all$YrSold <- as.factor(all$YrSold) #the numeric version is now not needed anymore


unnamed-chunk-105-1.png


7.3 Binning Neighborhood


R

nb1 <- ggplot(all[!is.na(all$SalePrice),], aes(x=reorder(Neighborhood, SalePrice, FUN=median), y=SalePrice)) +

geom_bar(stat='summary', fun.y = "median", fill='blue') + labs(x='Neighborhood', y='Median SalePrice') +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_y_continuous(breaks= seq(0, 800000, by=50000), labels = comma) +
geom_label(stat = "count", aes(label = ..count.., y = ..count..), size=3) +
geom_hline(yintercept=163000, linetype="dashed", color = "red") #dashed line is median SalePrice
nb2 <- ggplot(all[!is.na(all$SalePrice),], aes(x=reorder(Neighborhood, SalePrice, FUN=mean), y=SalePrice)) +
geom_bar(stat='summary', fun.y = "mean", fill='blue') + labs(x='Neighborhood', y="Mean SalePrice") +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_y_continuous(breaks= seq(0, 800000, by=50000), labels = comma) +
geom_label(stat = "count", aes(label = ..count.., y = ..count..), size=3) +
geom_hline(yintercept=163000, linetype="dashed", color = "red") #dashed line is median SalePrice
grid.arrange(nb1, nb2)

unnamed-chunk-107-1.png

中央値および平均Salepricesは、実質的に高い金額を持つ3つのNeighborhoodで一致します。3つの比較的貧しい地域の分離はあまり明確ではありませんが、少なくとも両方のグラフは同じ3つの貧しい地域が一致しています。


R

all$NeighRich[all$Neighborhood %in% c('StoneBr', 'NridgHt', 'NoRidge')] <- 2

all$NeighRich[!all$Neighborhood %in% c('MeadowV', 'IDOTRR', 'BrDale', 'StoneBr', 'NridgHt', 'NoRidge')] <- 1
all$NeighRich[all$Neighborhood %in% c('MeadowV', 'IDOTRR', 'BrDale')] <- 0

table(all$NeighRich)
##
## 0 1 2
## 160 2471 288



7.4 Total Square Feet

人々が住宅を購入するとき、全体の生活スペースは一般的に非常に重要なので、地上と地下の居住スペースを合わせた変数を追加しています。


R

all$TotalSqFeet <- all$GrLivArea + all$TotalBsmtSF

ggplot(data=all[!is.na(all$SalePrice),], aes(x=TotalSqFeet, y=SalePrice))+
geom_point(col='blue') + geom_smooth(method = "lm", se=FALSE, color="black", aes(group=1)) +
scale_y_continuous(breaks= seq(0, 800000, by=100000), labels = comma) +
geom_text_repel(aes(label = ifelse(all$GrLivArea[!is.na(all$SalePrice)]>4500, rownames(all), '')))

unnamed-chunk-111-1.png

予想どおり、SalePriceとの相関は確かに非常に強い(0.78)。


R

cor(all$SalePrice, all$TotalSqFeet, use= "pairwise.complete.obs")

## [1] 0.7789588

2つの潜在的な外れ値は、以前よりさらに「外れ値っぽい」と思われる。これら2つの外れ値を取り出すことによって、相関は5%増加します。


R

cor(all$SalePrice[-c(524, 1299)], all$TotalSqFeet[-c(524, 1299)], use= "pairwise.complete.obs")

## [1] 0.829042


7.5 Consolidating Porch variables

以下では、私はPorchに関連すると思われる変数を挙げました。


  • WoodDeckSF:木製デッキ面積(平方フィート)

  • OpenPorchSF:オープンポーチエリア(平方フィート)

  • EnclosedPorch:エンクローズポーチエリア(平方フィート)

  • 3SsnPorch:3シーズンのポーチエリア(平方フィート)

  • ScreenPorch:スクリーンポーチの面積(平方フィート)

私の知る限りでは、Porchは家の外で保護された場所であり、木製のデッキは隠されていません。したがって、私はWoodDeckSFだけを残しており、4つのポーチ変数を統合しています。


R

all$TotalPorchSF <- all$OpenPorchSF + all$EnclosedPorch + all$X3SsnPorch + all$ScreenPorch


これらのポーチのエリアを合計すると意味がありますが(領域間の重複はない)、SalePriceとの相関はそれほど強くはありません。


R

cor(all$SalePrice, all$TotalPorchSF, use= "pairwise.complete.obs")

## [1] 0.1957389
ggplot(data=all[!is.na(all$SalePrice),], aes(x=TotalPorchSF, y=SalePrice))+
geom_point(col='blue') + geom_smooth(method = "lm", se=FALSE, color="black", aes(group=1)) +
scale_y_continuous(breaks= seq(0, 800000, by=100000), labels = comma)

unnamed-chunk-116-1.png


8 Preparing data for modeling


8.1 Dropping highly correlated variables

まず、2つの変数の相関性が高い場合は変数を削除します。これらの相関ペアを見つけるために、相関行列を再度使用しました(6.1節を参照)。たとえば、GarageCarsとGarageAreaの相関は0.89です。これら2つのうち、SalePrice(GarageAreaとSalePrice相関が0.62、GarageCarsのSalePrice相関が0.64です)との相関が最も低い変数を削除します。


R

dropVars <- c('YearRemodAdd', 'GarageYrBlt', 'GarageArea', 'GarageCond', 'TotalBsmtSF', 'TotalRmsAbvGrd', 'BsmtFinSF1')

all <- all[,!(names(all) %in% dropVars)]



8.2 Removing outliers

暫くの間、外れ値を単純にキープしましたが、SalePriceが低い2つの外れ値を手動で取り除きます。しかし、後の段階でもっと徹底的に調べるつもりです。


R

all <- all[-c(524, 1299),]



8.3 PreProcessing predictor variables

モデリングの前に、数値の予測子(ラベルにエンコードされた変数ではない)をスケールし、カテゴリの予測変数のダミー変数を作成する必要があります。以下では、データフレームをすべて数値変数と、(序数)因子を保持する別のデータフレームに分割しています。


R

numericVarNames <- numericVarNames[!(numericVarNames %in% c('MSSubClass', 'MoSold', 'YrSold', 'SalePrice', 'OverallQual', 'OverallCond'))] #numericVarNames was created before having done anything

numericVarNames <- append(numericVarNames, c('Age', 'TotalPorchSF', 'TotBathrooms', 'TotalSqFeet'))

DFnumeric <- all[, names(all) %in% numericVarNames]

DFfactors <- all[, !(names(all) %in% numericVarNames)]
DFfactors <- DFfactors[, names(DFfactors) != 'SalePrice']

cat('There are', length(DFnumeric), 'numeric variables, and', length(DFfactors), 'factor variables')
## There are 30 numeric variables, and 49 factor variables



8.3.1 Skewness and normalizing of the numeric predictors

歪度は、分布における対称性の尺度です。対称なデータセットの歪度は0になります。したがって、正規分布の歪度は0になります。歪度は本質的に2つの裾の相対的なサイズを測定します。経験則として、歪度は-1と1の間でなければなりません。この範囲では、データは対称的と見なされます。歪度を修正するために、絶対歪度が0.8より大きいすべての数値予測子の対数を取っています(実際にはゼロ問題による除算を避けるためにlog + 1です)。


R

for(i in 1:ncol(DFnumeric)){

if (abs(skew(DFnumeric[,i]))>0.8){
DFnumeric[,i] <- log(DFnumeric[,i] +1)
}
}

PreNum <- preProcess(DFnumeric, method=c("center", "scale"))
print(PreNum)
## Created from 2917 samples and 30 variables
##
## Pre-processing:
## - centered (30)
## - ignored (0)
## - scaled (30)

DFnorm <- predict(PreNum, DFnumeric)
dim(DFnorm)
## [1] 2917 30



8.3.2 One hot encoding the categorical variables

すべての予測変数を数値列に変換するために必要な最後のステップは、カテゴリ変数を「ワンホットエンコード」することです。これは、変数を1と0で別々の列で表現することです。このOne hot encodingを行うには、model.matrix関数を使用しています。


R

DFdummies <- as.data.frame(model.matrix(~.-1, DFfactors))

dim(DFdummies)
## [1] 2917 201


8.3.3 Removing levels with few or no observations in train or test

以前は、CaretNear Zero Variance関数で作業しました。 これは機能しますが、多くの情報が失われました。たとえば、デフォルトを使用すると、146未満のすべてのNeighborhoodsが省略されます。 したがって、もっと慎重な手動なアプローチを取ってきました。


R

#check if some values are absent in the test set

ZerocolTest <- which(colSums(DFdummies[(nrow(all[!is.na(all$SalePrice),])+1):nrow(all),])==0)
colnames(DFdummies[ZerocolTest])
## [1] "Condition2RRAe" "Condition2RRAn" "Condition2RRNn"
## [4] "HouseStyle2.5Fin" "RoofMatlMembran" "RoofMatlMetal"
## [7] "RoofMatlRoll" "Exterior1stImStucc" "Exterior1stStone"
## [10] "Exterior2ndOther" "HeatingOthW" "ElectricalMix"
## [13] "MiscFeatureTenC"
DFdummies <- DFdummies[,-ZerocolTest] #removing predictors
#check if some values are absent in the train set
ZerocolTrain <- which(colSums(DFdummies[1:nrow(all[!is.na(all$SalePrice),]),])==0)
colnames(DFdummies[ZerocolTrain])
## [1] "MSSubClass1,5 story PUD all"
DFdummies <- DFdummies[,-ZerocolTrain] #removing predictor
Also taking out variables with less than 10 ones in the train set.

fewOnes <- which(colSums(DFdummies[1:nrow(all[!is.na(all$SalePrice),]),])<10)
colnames(DFdummies[fewOnes])
## [1] "MSSubClass1 story unf attic" "LotConfigFR3"
## [3] "NeighborhoodBlueste" "NeighborhoodNPkVill"
## [5] "Condition1PosA" "Condition1RRNe"
## [7] "Condition1RRNn" "Condition2Feedr"
## [9] "Condition2PosA" "Condition2PosN"
## [11] "RoofStyleMansard" "RoofStyleShed"
## [13] "RoofMatlWdShake" "RoofMatlWdShngl"
## [15] "Exterior1stAsphShn" "Exterior1stBrkComm"
## [17] "Exterior1stCBlock" "Exterior2ndAsphShn"
## [19] "Exterior2ndBrk Cmn" "Exterior2ndCBlock"
## [21] "Exterior2ndStone" "FoundationStone"
## [23] "FoundationWood" "HeatingGrav"
## [25] "HeatingWall" "ElectricalFuseP"
## [27] "GarageTypeCarPort" "MiscFeatureOthr"
## [29] "SaleTypeCon" "SaleTypeConLD"
## [31] "SaleTypeConLI" "SaleTypeConLw"
## [33] "SaleTypeCWD" "SaleTypeOth"
## [35] "SaleConditionAdjLand"
DFdummies <- DFdummies[,-fewOnes] #removing predictors
dim(DFdummies)
## [1] 2917 152


全体的には、分散がほとんどない49のワンホットコード予測子を削除しました。これはかなりの数に見えるかもしれませんが、実際にはcaretNear Zero Variance(デフォルトのしきい値を使用)を使用して取り出された予測変数の数よりもはるかに少ないです。


R

combined <- cbind(DFnorm, DFdummies) #combining all (now numeric) predictors into one dataframe 



8.4 Dealing with skewness of response variable


R

skew(all$SalePrice)

## [1] 1.877427
qqnorm(all$SalePrice)
qqline(all$SalePrice)

unnamed-chunk-129-1.png

歪度1.87の右の歪みが高すぎることを示し、QQプロットは販売価格もまた正規分布していないことを示しています。これを修正するため、私はSalePriceの対数を取っています。


R

all$SalePrice <- log(all$SalePrice) #default is the natural logarithm, "+1" is not necessary as there are no 0's

skew(all$SalePrice)
## [1] 0.1213182
As you can see,the skew is now quite low and the Q-Q plot is also looking much better.

qqnorm(all$SalePrice)
qqline(all$SalePrice)


unnamed-chunk-131-1.png


8.5 Composing train and test sets


R

train1 <- combined[!is.na(all$SalePrice),]

test1 <- combined[is.na(all$SalePrice),]


9 Modeling


9.1 Lasso regression model

私はRidgeとElastic Netのモデルも試しましたが、Lasso回帰がこれらの3つのモデルの中で最良の結果を出したので、ドキュメントにLasso回帰モデルだけを残しています。

Elastic Netのペナルティはαによって制御され、Lasso回帰(α= 1)とRidge回帰(α=0)との間の橋渡しをします。チューニングパラメータλは、ペナルティの全体的な強さを制御します。Ridge回帰は予測子の係数を互いに向かって縮小するが、一方、Lasso回帰はそれらの1つを選択して他のものを捨てる傾向があることが知られている。

以下では、Lasso回帰モデルに対してチューニングする必要がある唯一のハイパーパラメータであるλの最適値を見つけるために、caretのクロスバリデーションを使用しています。


R

set.seed(27042018)

my_control <-trainControl(method="cv", number=5)
lassoGrid <- expand.grid(alpha = 1, lambda = seq(0.001,0.1,by = 0.0005))

lasso_mod <- train(x=train1, y=all$SalePrice[!is.na(all$SalePrice)], method='glmnet', trControl= my_control, tuneGrid=lassoGrid)
lasso_mod$bestTune
## alpha lambda
## 4 1 0.0025
min(lasso_mod$results$RMSE)
## [1] 0.1121579


caretvarImp関数のドキュメントによると、glmboostglmnetの場合、チューニングされたモデルに対応する係数の絶対値が使用されます。

これは、最も重要な変数の実際のランクが格納されていないことを意味しますが、モデルで使用されていない変数の数(したがって係数は0)を調べる機会を与えてくれます。


R

lassoVarImp <- varImp(lasso_mod,scale=F)

lassoImportance <- lassoVarImp$importance

varsSelected <- length(which(lassoImportance$Overall!=0))
varsNotSelected <- length(which(lassoImportance$Overall==0))

cat('Lasso uses', varsSelected, 'variables in its model, and did not select', varsNotSelected, 'variables.')
## Lasso uses 100 variables in its model, and did not select 82 variables.


Lasso回帰はモデル内の利用可能な変数の約45%を使用しないことで、マルチコをうまく処理したようです。


R

LassoPred <- predict(lasso_mod, test1)

predictions_lasso <- exp(LassoPred) #need to reverse the log to the real values
head(predictions_lasso)
## 1461 1462 1463 1464 1465 1466
## 114351.8 162204.8 179455.3 197564.7 205952.8 169839.8


9.2 XGBoost model

最初は、XGBoostパッケージを使って直接作業します。理由は、パッケージが独自の効率的なデータ構造(xgb.DMatrix)を使用しているためです。パッケージには、CVの機能もあります。ただし、このCV関数は最適なラウンド数を決定するだけで、ハイパーパラメータのフルグリッドサーチはサポートしません。

caretxgbパッケージの(高速な)データ構造を使用しているようには見えませんが、少なくともフル・グリッド・サーチをサポートしているので、caretでハイパーパラメータ・チューニングを行うことにしました。私が理解する限り、オーバーフィットを避けるために調整する主なパラメータはmax_depthmin_child_weightです(XGBoostのドキュメントを参照)。 以下では、これらのパラメータとeta(学習率)の両方を調整するグリッドを設定しています。


R

xgb_grid = expand.grid(

nrounds = 1000,
eta = c(0.1, 0.05, 0.01),
max_depth = c(2, 3, 4, 5, 6),
gamma = 0,
colsample_bytree=1,
min_child_weight=c(1, 2, 3, 4 ,5),
subsample=1
)

次のステップは、caretで(5holdsCVを使用)最良のハイパーパラメータ値を見つけることです。


R

#xgb_caret <- train(x=train1, y=all$SalePrice[!is.na(all$SalePrice)], method='xgbTree', trControl= my_control, tuneGrid=xgb_grid) 

#xgb_caret$bestTune

予想どおり、これはかなり時間がかかりました。Kaggleの実行時間を制限したいので、コードを無効にして、結果を続行しています。caretによると、bestTuneのパラメータは次のとおりです。


  • Max_depth = 3

  • η = 0.05

  • Min_child_weight = 4

このセクションの残りの部分では、私はxgboostパッケージを使って直接作業を続けます。


R

label_train <- all$SalePrice[!is.na(all$SalePrice)]

# put our testing & training data into two seperates Dmatrixs objects
dtrain <- xgb.DMatrix(data = as.matrix(train1), label= label_train)
dtest <- xgb.DMatrix(data = as.matrix(test1))

default_param<-list(
objective = "reg:linear",
booster = "gbtree",
eta=0.05, #default = 0.3
gamma=0,
max_depth=3, #default=6
min_child_weight=4, #default=1
subsample=1,
colsample_bytree=1
)


のステップは、クロスバリデーションを実行して、(パラメータのセットに対して)最良のラウンド数を決定することです。


R

xgbcv <- xgb.cv( params = default_param, data = dtrain, nrounds = 500, nfold = 5, showsd = T, stratified = T, print_every_n = 40, early_stopping_rounds = 10, maximize = F)

## [1] train-rmse:10.955588+0.004477 test-rmse:10.955537+0.019106
## Multiple eval metrics are present. Will use test_rmse for early stopping.
## Will train until test_rmse hasn't improved in 10 rounds.
##
## [41] train-rmse:1.428274+0.000561 test-rmse:1.428515+0.011791
## [81] train-rmse:0.219833+0.000801 test-rmse:0.230612+0.009490
## [121] train-rmse:0.102497+0.001285 test-rmse:0.128856+0.008654
## [161] train-rmse:0.090461+0.001218 test-rmse:0.122128+0.007505
## [201] train-rmse:0.084142+0.001214 test-rmse:0.119557+0.007344
## [241] train-rmse:0.079398+0.001195 test-rmse:0.118374+0.007088
## [281] train-rmse:0.075716+0.001302 test-rmse:0.117645+0.006772
## [321] train-rmse:0.072567+0.001139 test-rmse:0.117136+0.006720
## [361] train-rmse:0.069770+0.001156 test-rmse:0.116745+0.006603
## [401] train-rmse:0.067201+0.001059 test-rmse:0.116505+0.006574
## [441] train-rmse:0.064814+0.001145 test-rmse:0.116386+0.006366
## Stopping. Best iteration:
## [454] train-rmse:0.063958+0.001067 test-rmse:0.116289+0.006326

これは少しの時間の作業でしたが、クロスバリデーションされたRMSEがかなり改善されているため、ハイパーパラメータチューニングが確実に行われています。


R

#train the model using the best iteration found by cross validation

xgb_mod <- xgb.train(data = dtrain, params=default_param, nrounds = 454)
XGBpred <- predict(xgb_mod, dtest)
predictions_XGB <- exp(XGBpred) #need to reverse the log to the real values
head(predictions_XGB)
## [1] 116386.8 162307.3 186494.0 187440.4 187258.3 166241.4
#view variable importance plot
library(Ckmeans.1d.dp) #required for ggplot clustering
mat <- xgb.importance (feature_names = colnames(train1),model = xgb_mod)
xgb.ggplot.importance(importance_matrix = mat[1:20], rel_to_first = TRUE)

unnamed-chunk-143-1.png


9.3 Averaging predictions

Lasso回帰とXGBoostのアルゴリズムは非常に異なるため、予測結果の平均化ではスコアが向上する可能性があります。Lasso回帰モデルがクロスバリデーションされたRMSEスコア(0.1121対0.1162)に関してより良くなるので、私はLasso回帰モデルに2倍を重み付けしています。


R

sub_avg <- data.frame(Id = test_labels, SalePrice = (predictions_XGB+2*predictions_lasso)/3)

head(sub_avg)
## Id SalePrice
## 1461 1461 115030.1
## 1462 1462 162238.9
## 1463 1463 181801.5
## 1464 1464 194189.9
## 1465 1465 199721.3
## 1466 1466 168640.3