R
RStudio
Kaggle

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