Help us understand the problem. What is going on with this article?

# kaggleのkernels学習ノート~Extensive Sberbank Exploratory Analysis~

More than 1 year has passed since last update.

# Extensive Sberbank Exploratory Analysis

• Missing Data
• Data Quality Issues
• Internal Home Features
• Demographic Features
• Education Features
• Cultural/Recreational Features
• Infrastructure Features
• Variable Importance

• Missing Data
• Internal Features
R
```library(data.table)
library(tidyverse)
library(lubridate)
library(scales)
library(corrplot)
library(DT)
```

## 訓練データ

R
```dtrain <- fread('../input/train.csv', stringsAsFactors=TRUE)
```

### 欠損データ

どのくらい欠損値があるのでしょうか。

R
```miss_pct <- map_dbl(dtrain, function(x) { round((sum(is.na(x)) / length(x)) * 100, 1) })

miss_pct <- miss_pct[miss_pct > 0]

data.frame(miss=miss_pct, var=names(miss_pct), row.names=NULL) %>%
ggplot(aes(x=reorder(var, -miss), y=miss)) +
geom_bar(stat='identity', fill='red') +
labs(x='', y='% missing', title='Percent missing data by feature') +
theme(axis.text.x=element_text(angle=90, hjust=1))
```

292列のうち、51列は欠損値を持っています。欠落率は、metro_min_walkの0.1％からhospital_beds_raionで47.4％までの範囲です。

### データ品質の問題

R
```# state should be discrete valued between 1 and 4. There is a 33 in it that is cleary a data entry error
# Lets just replace it with the mode.
dtrain\$state[dtrain\$state == 33] <- which.max(table(dtrain\$state))

# build_year has an erronus value 20052009. Since its unclear which it should be, let's replace with 2007
dtrain\$build_year[dtrain\$build_year == 20052009] <- 2007
```

### 家の内部の特徴

```R:R# variables that are internal to home
internal_chars <- c('full_sq', 'life_sq', 'floor', 'max_floor', 'build_year', 'num_room',
'kitch_sq', 'state', 'price_doc')

corrplot(cor(dtrain[, ..internal_chars], use="complete.obs"))
```

### 家の面積と部屋の数

R
```ggplot(aes(x=full_sq, y=price_doc), data=dtrain) +
geom_point(color='red')
```

`full_sq`に異常値があります。これがエントリーエラーかどうかは不明です。今のところは除外します。

R
```dtrain %>%
filter(full_sq < 2000) %>%
ggplot(aes(x=full_sq, y=price_doc)) +
geom_point(color='red', alpha=0.5) +
labs(x='Area', y='Price', title='Price by area in sq meters')
```

`full_sq`は、データ・ディクショナリーでは、バルコニー、その他の非居住地域を含む総面積（平方メートル）で定義され、`life_sq`は、「バルコニー、その他の非居住地域を除く平方メートル単位の生活面積」として定義されています。したがって、`life_sq`は常に`full_sq`未満である必要があります。

R
```sum(dtrain\$life_sq > dtrain\$full_sq, na.rm=TRUE)
## [1] 37
```

`life_sq``full_sq`より大きい観測値が37個あります。`num_rooms`を見てみましょう。

R
```table(dtrain\$num_room)
##
##    0    1    2    3    4    5    6    7    8    9   10   17   19
##   14 7602 8132 4675  418   40    9    1    3    1    2    1    1
ggplot(aes(x=num_room), data=dtrain) +
geom_histogram(fill='red', bins=20) +
ggtitle('Distribution of room count')
## Warning: Removed 9572 rows containing non-finite values (stat_bin).
```

### 販売タイプ

R
```ggplot(aes(x=price_doc), data=dtrain) +
geom_density(fill='red', color='red') +
facet_grid(~product_type) +
scale_x_continuous(trans='log')
```

R
```dtrain %>% group_by(product_type) %>% summarize(median(price_doc))
## # A tibble: 2 × 2
##    product_type `median(price_doc)`
##          <fctr>               <dbl>
## 1    Investment             6670000
## 2 OwnerOccupier             5564090
```

### 建築年

`build_year`を見てみましょう。

R
```table(dtrain\$build_year)
##
##    0    1    3   20   71  215 1691 1860 1876 1886 1890 1895 1896 1900 1904
##  530  368    2    1    1    1    1    2    1    1    5    1    2    2    1
## 1905 1906 1907 1910 1911 1912 1914 1915 1917 1920 1924 1925 1926 1927 1928
##    1    1    2    5    1    5    3    5   16    1    3    1    8   10   12
## 1929 1930 1931 1932 1933 1934 1935 1936 1937 1938 1939 1940 1941 1943 1946
##   12    6    6    8    7   13   11    5   10    9    7   14    2    2    2
## 1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961
##    4    1    3   22   23   45   24   36   52   48  119  179  208  344  297
## 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976
##  338  325  315  378  348  384  389  407  418  352  360  333  357  309  263
## 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991
##  260  235  236  226  189  189  185  169  178  131  171  155  155  129   93
## 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006
##  139  115  160  149  162  139  141  125  130  177  214  193  220  176  242
## 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 4965
##  220  234  176  132  162  233  464  919  824  375  154    1    1

```

R
```dtrain %>%
filter(build_year > 1691 & build_year < 2018) %>%
ggplot(aes(x=build_year)) +
geom_histogram(fill='red') +
ggtitle('Distribution of build year')
```

この分布は、1970年代初期と過去数年間にピークを持つ二峰性のように見えます。今、`build_year``price`が関連しているかどうかを見てみましょう。 ここでは、年ごとにデータをグループ化し、`price_doc`の平均値を取っています。

R
```dtrain %>%
filter(build_year > 1691 & build_year < 2018) %>%
group_by(build_year) %>%
summarize(mean_build_price=mean(price_doc)) %>%
ggplot(aes(x=build_year, y=mean_build_price)) +
geom_line(stat='identity', color='red') +
geom_smooth(color='darkgrey') +
ggtitle('Mean price by year of build')
```

この関係は、特に1960年以降、安定しているように見えます。初期の年代にはいくつかの変動があります。これは実際の効果ではなく、単に1950年頃までの観察の希薄さによるものです。

### タイムスタンプ

R
```dtrain\$timestamp <- as.Date(dtrain\$timestamp)

dtrain %>%
group_by(timestamp) %>%
summarize(med_price = median(price_doc)) %>%
ggplot(aes(x = timestamp, y = med_price)) +
geom_line(color = 'red') +
geom_smooth(method = 'lm', color = 'grey', alpha = 0.7) +
ggtitle('Daily median price over time')
```

そして上記のプロットと比較すると、同じ時間に売り上げのボリュームがあります。

R
```dtrain %>%
group_by(timestamp) %>%
summarize(n = n()) %>%
ggplot(aes(x = timestamp, y = n)) +
geom_bar(stat = 'identity') +
labs(x='', y='Number of transactions', title='Sales volume over time')
```

1年の間に住宅価格に季節要因がありますか？

R
```dtrain %>%
mutate(month=month(timestamp)) %>%
group_by(month) %>%
summarize(med_price=median(price_doc)) %>%
ggplot(aes(x=as.integer(month), y=med_price)) +
geom_line(color='red', stat='identity') +
geom_point(color='red', size=2) +
scale_x_continuous(breaks=seq(1,12,1)) +
labs(x='Month', title='Price by month of year')
```

### 家の状態や材質

R
```dtrain %>%
filter(!is.na(state)) %>%
ggplot(aes(x=as.factor(state), y=log10(price_doc))) +
geom_jitter(color='grey', alpha=0.2) +
geom_violin(fill='red', alpha=0.7) +
ggtitle('Log10 of median price by state of home')
```

プロットからでは分かりにくいですが、`state==4`は平均で最も高い販売価格を持っているようです。しかし、このカテゴリーは、家の数がかなり少ないです。確認しよう。

R
```dtrain %>%
filter(!is.na(state)) %>%
group_by(state) %>%
summarize(mean(price_doc))
## # A tibble: 4 × 2
##   state `mean(price_doc)`
##   <int>             <dbl>
## 1     1           7315440
## 2     2           7060396
## 3     3           8078316
## 4     4          13345469
```

`state==4`は平均価格が最も高く、`state==3`に続きます。`state==1``state==2`は近い状態です。材質の特徴量はどうでしょうか？

R
```table(dtrain\$material)
##
##     1     2     3     4     5     6
## 14197  2993     1  1344  1561   803
```

この特徴量はデータディクショナリーには記述されていないため、これらの値が何を意味するのかは不明です。`material==1`は一般的です。 1つの家だけが`material==3`として分類されます。これら6つの材料の中央値はどのようにして比較されますか？

R
```dtrain %>%
filter(!is.na(material)) %>%
ggplot(aes(x=as.factor(material), y=log(price_doc))) +
geom_jitter(alpha=0.4, color='grey') +
geom_violin(fill='red', color='red',  alpha=0.6) +
ggtitle('Distribution of price by build material')
```

R
```dtrain %>%
filter(!is.na(material)) %>%
group_by(state=as.factor(material)) %>%
summarize(med_price=median(price_doc))
## # A tibble: 6 × 2
##    state med_price
##   <fctr>     <dbl>
## 1      1   6500000
## 2      2   6900000
## 3      3   6931143
## 4      4   7247870
## 5      5   6492000
## 6      6   6362318
```

それほど大きな違いはありませんが、`state==4`は最も高価なようです。

### 家のフロア

フロアは価格を比べるとどうでしょうか。以前の相関プロットによると、中程度の正の相関があります。

R
```ggplot(aes(x=floor, y=log(price_doc)), data=dtrain) +
geom_point(color='red', alpha=0.4) +
geom_smooth(method='lm', color='darkgrey') +
ggtitle('Price by floor of home')
```

R
```ggplot(aes(x=max_floor, y=log(price_doc)), data=dtrain) +
geom_point(color='red', alpha=0.4) +
geom_smooth(method='lm', color='darkgrey')
```

データの品質をチェックし、床がどの観測値でも`max_floor`より大きいかどうかを見てみましょう。

R
```ggplot(aes(x=floor, y=max_floor), data=dtrain) +
geom_point(color='red') +
geom_abline(slope=1, intercept=0, color='grey')
## Warning: Removed 9572 rows containing missing values (geom_point).
```

R
```dtrain %>%
select(id, floor, max_floor) %>%
filter(floor > max_floor) %>%
datatable()
```
id floor max_floor
1 8219 13 0
2 8271 3 1
3 8502 2 0
4 8534 7 0
5 8915 5 0
6 9164 8 3
7 9260 8 1
8 9312 5 1
9 9391 10 1
10 9415 4 1

id floor max_floor
1491 30429 1 0
1492 30442 12 0
1493 30453 5 0

この場合、1,493件の観測があります。

## デモグラフィック情報

ここでは、家の内部の特性を超えて、基本的なデモグラフィックとジオグラフィックを見てみましょう。まず、相関プロットです。

R
```demo_vars <- c('area_m', 'raion_popul', 'full_all', 'male_f', 'female_f', 'young_all',
'young_female', 'work_all', 'work_male', 'work_female', 'price_doc')

corrplot(cor(dtrain[, ..demo_vars], use='complete.obs'))
```

R
```# How many unique districts are there?
length(unique(dtrain\$sub_area))
## [1] 146
```

R
```dtrain %>%
mutate(area_km=area_m/1000000, density=raion_popul/area_km) %>%
select(sub_area, density, price_doc) %>%
group_by(sub_area) %>%
summarize(density=median(density), med_price=median(price_doc)) %>%
ggplot(aes(x=density, y=med_price)) +
geom_point(color='grey') +
geom_smooth(method='lm', color='red') +
ggtitle('Median home price by raion population density (people per sq. km)')
```

これらの数字は、モスクワ全体の人口密度が8,537/sqkmであることを考えると意味をなさないようです。奇妙な値に見えるいくつかの地区があります。 住宅価格は人口密度とともに増加するようです。各地区に販売取引がいくつあるのか考えてみましょう。

R
```dtrain %>%
group_by(sub_area) %>%
summarize(n=n()) %>%
ggplot(aes(x=reorder(sub_area, n), y=n)) +
geom_bar(stat='identity') +
coord_flip() +
labs(y='Number of transactions', x='', title='Number of Transactions by District')
```

Poselenie Sosenskoe、Nekrasovka、Poselenie Vnukovskoeは、データセット内で最も多く取引を行っていました。働く年齢と価格の人口の割合は、関係があると思います。

R
```dtrain %>%
mutate(work_share=work_all/raion_popul) %>%
group_by(sub_area) %>%
summarize(mean_price=mean(price_doc), work_share=mean(work_share)) %>%
ggplot(aes(x=work_share, y=mean_price)) +
geom_point(color='red') +
geom_smooth(color='gray') +
ggtitle('District mean home price by share of working age population')
```

### School Characteristics

R
```school_chars <- c('children_preschool', 'preschool_quota', 'preschool_education_centers_raion',
'children_school', 'school_quota', 'school_education_centers_raion',
'school_education_centers_top_20_raion', 'university_top_20_raion',
'price_doc')

corrplot(cor(dtrain[, ..school_chars], use='complete.obs'))
```

R
```table(dtrain\$university_top_20_raion)
##
##     0     1     2     3
## 27392  1994  1035    50
dtrain %>%
ggplot(aes(x=as.factor(university_top_20_raion), y=price_doc)) +
geom_jitter(color='grey') +
geom_boxplot(fill='red', color='gray', alpha=0.5) +
ggtitle('Distribution of home price by # of top universities in Raion')
```

3つのトップ20の大学を持つ地区の家庭は、最高の住宅価格の中央値を持っていますが、0,1,2の間でかなり近いです。3つのトップ大学を持つ家がほとんどありません。 3つの大学といくつの地区があるのか見てみましょう。

R
```unique(dtrain %>% filter(university_top_20_raion==3) %>% select(sub_area))
##         sub_area
## 1 Zamoskvorech'e
```

## 文化/レクリエーションの特性

R
```cult_chars <- c('sport_objects_raion', 'culture_objects_top_25_raion', 'shopping_centers_raion',                        'park_km', 'fitness_km', 'swim_pool_km', 'ice_rink_km','stadium_km', 'basketball_km',                   'shopping_centers_km', 'big_church_km','church_synagogue_km', 'mosque_km', 'theater_km',                 'museum_km', 'exhibition_km', 'catering_km', 'price_doc')

corrplot(cor(dtrain[, ..cult_chars], use='complete.obs'))
```

R
```dtrain %>%
group_by(sub_area) %>%
summarize(sport_objects=mean(sport_objects_raion), med_price=median(price_doc)) %>%
ggplot(aes(x=sport_objects, y=med_price)) +
geom_point(color='grey') +
geom_smooth(method='lm', color='red') +
ggtitle('Median Raion home price by # of sports objects in Raion')
```

R
```dtrain %>%
group_by(sub_area) %>%
summarize(culture_objects=mean(culture_objects_top_25_raion), med_price=median(price_doc)) %>%
ggplot(aes(x=culture_objects, y=med_price)) +
geom_point(color='grey') +
geom_smooth(method='lm', color='red') +
ggtitle('Median raion home price by # of culture objects in raion')
```

R
```dtrain %>% group_by(culture_objects_top_25) %>%
summarize(med_price=median(price_doc))
## # A tibble: 2 × 2
##   culture_objects_top_25 med_price
##                   <fctr>     <dbl>
## 1                     no   6200000
## 2                    yes   7400000

```

したがって、上位25の文化的対象を持つ地区は、120万人以上の販売価格の中央値を持っています。最も近い公園までの距離は住宅価格にどのように関係していますか？

R
```ggplot(aes(x=park_km, y=price_doc), data=dtrain) +
geom_point(color='red', alpha=0.4) +
geom_smooth(method='lm', color='grey') +
ggtitle('Home price by distance to nearest park')
```

## インフラ特性

R
```inf_features <- c('nuclear_reactor_km', 'thermal_power_plant_km', 'power_transmission_line_km',

corrplot(cor(dtrain[, ..inf_features], use='complete.obs'))
```

R
```ggplot(aes(x=kremlin_km, y=price_doc), data=dtrain) +
geom_point(color='grey') +
geom_smooth(method='lm', color='red') +
ggtitle('Home price by distance to Kremlin')
```

## 変数重要度

R
```library(caret)

#Get complete cases of dtrain
completes <- complete.cases(dtrain)

# Set training control so that we only 1 run forest on the entire set of complete cases
trControl <- trainControl(method='none')

# Run random forest on complete cases of dtrain. Exclude incineration_raion since it
# only has 1 factor level
rfmod <- train(price_doc ~ . - id - timestamp - incineration_raion,
method='rf',
data=dtrain[completes, ],
trControl=trControl,
tuneLength=1,
importance=TRUE)

varImp(rfmod)
## rf variable importance
##
##   only 20 most important variables shown (out of 435)
##
##                             Overall
## full_sq                      100.00
## life_sq                       86.79
## num_room                      86.77
## kitch_sq                      74.05
## build_year                    67.96
## max_floor                     54.31
## product_typeOwnerOccupier     51.06
## state                         40.94
## theater_km                    38.91
## metro_km_avto                 38.63
## workplaces_km                 37.92
## metro_min_avto                37.80
## cafe_sum_3000_min_price_avg   37.75
## cemetery_km                   37.62
## hospice_morgue_km             37.59
## sport_count_2000              37.57
## ice_rink_km                   37.53
```

## Train vs Test Data

フォーラムでは、このコンテストのテストデータに関するいくつかの特徴量のディストリビューションが大きく異なっています。これにより、ローカルのクロスバリデーションスコアとリーダーボードスコアとの間に大きな相違が生じることがあります。 私たちが見つけることができるものを見てみましょう。

R
```dtest <- fread('../input/test.csv', stringsAsFactors=TRUE)
```

### 欠損値

R
```miss_pct <- map_dbl(dtest, function(x) { round((sum(is.na(x)) / length(x)) * 100, 1) })

miss_pct <- miss_pct[miss_pct > 0]

data.frame(miss=miss_pct, var=names(miss_pct), row.names=NULL) %>%
ggplot(aes(x=reorder(var, -miss), y=miss)) +
geom_bar(stat='identity', fill='red') +
labs(x='', y='% missing', title='Percent missing data by feature') +
theme(axis.text.x=element_text(angle=90, hjust=1))
```

R
```dtrain <- dtrain %>%
# remove price_doc from dtrain
select(-price_doc) %>%
mutate(dataset='train')

dtest <- dtest %>%
mutate(dataset='test', timestamp=as.Date(timestamp))

all_data <- bind_rows(dtrain, dtest)
```

### 室内の特徴

R
```all_data %>%
ggplot(aes(x=full_sq)) +
geom_density(color='red', fill='red', alpha=0.7) +
facet_wrap(~as.factor(dataset)) +
scale_x_continuous(trans='log') +
ggtitle('Distribution of full_sq')
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Removed 3 rows containing non-finite values (stat_density).
```

R
```all_data %>%
ggplot(aes(x=life_sq)) +
geom_density(color='red', fill='red', alpha=0.7) +
facet_wrap(~as.factor(dataset)) +
scale_x_continuous(trans='log') +
ggtitle('Distribution of life_sq')
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Removed 7608 rows containing non-finite values (stat_density).
```

R
```all_data %>%
ggplot(aes(x=kitch_sq)) +
geom_density(color='red', fill='red', alpha=0.7) +
facet_wrap(~as.factor(dataset)) +
scale_x_continuous(trans='log') +
ggtitle('Distribution of kitch_sq')
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Removed 11329 rows containing non-finite values (stat_density).
```

3つのエリアの分布は、訓練データとテストセットとの間で同様に分布していると思われますが、部屋数はどうでしょうか？

R
```all_data %>%
ggplot(aes(x=num_room)) +
geom_histogram(fill='red') +
facet_wrap(~dataset) +
ggtitle('Distribution of number of rooms')
## Warning: Removed 9572 rows containing non-finite values (stat_bin).
```

R
```all_data %>%
ggplot(aes(x=floor)) +
geom_density(color='red', fill='red') +
facet_wrap(~dataset)
## Warning: Removed 167 rows containing non-finite values (stat_density).
```

R
```all_data %>%
ggplot(aes(x=max_floor)) +
geom_density(color='red', fill='red') +
facet_wrap(~dataset)
## Warning: Removed 9572 rows containing non-finite values (stat_density).
```

いいですね。ここで`floor``max_floor`のデータを見て、`floor``max_floor`より大きいかどうかを確認します。

R
```all_data %>%
ggplot(aes(x=floor, y=max_floor)) +
geom_point(color='red') +
geom_abline(intercept=0, slope=1, color='darkgrey') +
facet_wrap(~dataset) +
ggtitle('Max floor by floor')
## Warning: Removed 9572 rows containing missing values (geom_point).
```

データの構造は訓練データとテストデータ間で同じ比率で持つようです。誤ったフロアデータの構造をより深く見て、構造があるかどうかを調べることが重要です。次に、訓練データとテストデータでトランザクションの数がどのように変化するかを見てみましょう。

R
```all_data %>%
ggplot(aes(x=timestamp, fill=dataset, color=dataset)) +
geom_bar(alpha=0.7) +
scale_fill_manual(values=c('red', 'darkgrey')) +
scale_color_manual(values=c('red', 'darkgrey')) +
ggtitle('Number of transactions by day')
```

R
```ggplot(aes(x=product_type),data=all_data) +
geom_bar(fill='red') +
facet_grid(~dataset)
```

テストセットには、`product_type`の値が欠落しています。最も一般的なタイプであるため、これらを「Investment」に置き換えることができます。販売数は、訓練データとテストデータ間で家の状態によってどのように異なるのでしょうか。

R
```all_data %>%
ggplot(aes(x=as.factor(state), fill=as.factor(state))) +
geom_bar() +
facet_wrap(~dataset) +
ggtitle('Distribution of state')
```

R
```all_data %>%
ggplot(aes(x=as.factor(material), fill=as.factor(material))) +
geom_bar() +
facet_wrap(~dataset) +
ggtitle('Distribution of material')
```

Why not register and get more from Qiita?
1. We will deliver articles that match you
By following users and tags, you can catch up information on technical fields that you are interested in as a whole
2. you can read useful information later efficiently
By "stocking" the articles you like, you can search right away