Функции и модули

mymean <- function(x) {
  return(sum(x) / length(x))
}
mysd <- function(x) {
  return(sqrt(sum((mymean(x) - x)^2)/(length(x) - 1)))
}
se_two <- function(x, y, type="general", sx=1, sy=2) {
  nx = length(x)
  ny = length(y)
  if (type == "general") {
    return(sqrt(sx^2/nx + sy^2 / ny))
  }
  if (type == "equal") {
    return(sx*sqrt(1 / nx + 1 / ny))
  }
  if (type == "unknown") {
    s = sqrt(((nx - 1) * mysd(x)^2 + (ny - 1) * mysd(y)^2)/ (nx + ny - 2))
    return(s * sqrt(1 / nx + 1 / ny))
  }
}
my_df_t <- function(x, y, type="general") {
  nx = length(x)
  ny = length(y)
  if (type == "general") {
    sx = mysd(x)
    sy = mysd(y)
    return((sx^2/nx + sy^2/ny)^2/((sx^2/nx)^2/(nx - 1) + (sy^2/ny)^2/(ny - 1)))
  }
  if (type == "equal") {
    return(nx + ny - 2)
  }
  if (type == "paired") {
    return(nx - 1)
  }
}
transform_func <- function(x, value) {
  if (x < value) {
    return(paste("<", value, sep=" "))
  } 
  else {
    return(paste(">=", value, sep=" "))
  }
}
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(reshape2)
library(RColorBrewer)
library(caTools)

Подготовим набор данных.

houses <- read.csv("housing.csv", header = TRUE, sep=",")
houses <- na.omit(houses)
unique(houses$ocean_proximity)
## [1] NEAR BAY   <1H OCEAN  INLAND     NEAR OCEAN ISLAND    
## Levels: <1H OCEAN INLAND ISLAND NEAR BAY NEAR OCEAN
houses <- mutate(houses, log10_value = log10(median_house_value))
house_list <- split(houses, houses$ocean_proximity)
size <- dim(houses)[1]
mean_value <- sum(houses$log10_value) / size
sd_value = sqrt(sum((mean_value - houses$log10_value)^2)/(size - 1))

Сравнение цен на нормальность

Предполагая стоимость домов распределённой логнормально, посмотрим на qqplot:

set.seed(123)
normal_q <- rnorm(size, mean = mean_value, sd = sd_value)
qqnorm(houses$log10_value)
qqline(normal_q)

Посмотрим на гистограмму:

hist(log10(houses$median_house_value), prob=TRUE)
curve(dnorm(x, mean=mean_value, sd = sd_value), xlim=c(4,6), add=TRUE, col='red')

Тест Колмогорова-Смирнова:

ks.test(x = log10(houses$median_house_value), y = normal_q)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  log10(houses$median_house_value) and normal_q
## D = 0.03416, p-value = 8.825e-11
## alternative hypothesis: two-sided

p-value мало, а значит, распределение не является нормальным. По графику видим, что есть небольшое утяжеление правого хвоста, что, по-видимому, обусловлено распространённостью элитной недвижимости.

Тест 1

H0: Нет значимой разницы между стоимостями домов на материке и рядом с бухтой. 

H1: Есть значимая разница на уровне 0.05. 

Логарифм стоимости распределён почти нормально, большая выборка, ничего не знаем про среднее и дисперсию. 

Нормальный тест.

z_val = (mymean(house_list$`NEAR BAY`$log10_value) - mymean(house_list$INLAND$log10_value)) / se_two(house_list$`NEAR BAY`$log10_value, house_list$INLAND$log10_value, type="unknown")
p_val <- pnorm(abs(z_val), lower.tail = FALSE) + pnorm(-abs(z_val))
p_val
## [1] 7.097714e-177

p-value < 0.05, поэтому отвергаем гипотезу H0.

Тест 2

H0: Нет значимой разницы между стоимостями домов на острове и рядом с бухтой.  

H1: стоимость домов на острове больше стоимости домов рядом с бухтой на уровне значимости 0.05.

qqplot(house_list$`NEAR BAY`$log10_value, house_list$ISLAND$log10_value)

t_val <- (mymean(house_list$ISLAND$log10_value) - mymean(house_list$`NEAR BAY`$log10_value)) / se_two(house_list$ISLAND$log10_value, house_list$`NEAR BAY`$log10_value, type="unknown")
df_val <- my_df_t(house_list$ISLAND$log10_value, house_list$`NEAR BAY`$log10_value, type="general")

Рассматриваем правостороннюю гипотезу

p_val <- pt(t_val, df=df_val, lower.tail = FALSE)
p_val
## [1] 0.1858747

Не отвергаем гипотезу H0.  

Посмотрим на нормальный тест.

z_val = (mymean(house_list$ISLAND$log10_value) - mymean(house_list$`NEAR BAY`$log10_value)) / se_two(house_list$ISLAND$log10_value, house_list$`NEAR BAY`$log10_value, type="unknown")
p_val <- pnorm(abs(z_val), lower.tail = FALSE)
p_val
## [1] 0.1581158

Вновь не отвергаем гипотезу H0.

Тест 3

H0: дома распределены равномерно по дальности от океана.  

H1: дома значимо распределены по-разному от океана на уровне значимости 0.01.

n_class <- length(unique(houses$ocean_proximity))
df <- n_class - 1
total <- dim(houses)[1]
expected <- total / n_class
observed <- as.vector(table(houses$ocean_proximity))
chi_stat <- sum((observed - expected)^2 / expected)
p_val <- pchisq(chi_stat, df = df, lower.tail = FALSE)
p_val
## [1] 0

На уровне значимости 0.01 отвергаем H0.

Тест 4

hist(houses$latitude)

hist(houses$longitude)

H0: широта и долгота домов независимы В группах “широта < 36” и “широта >=36” и “долгота < -120” и “долгота >= -120”.  

H1: широта и долгота домов зависимы в этих группах.

longitude_factor <- sapply(X = houses$longitude, FUN=transform_func, value=-120)
latitiude_factor <- sapply(X = houses$latitude, FUN=transform_func, value=36)
df_chisq <- (2 - 1) * (2 - 1)
observed <- as.matrix(table(data.frame(longitude_factor, latitiude_factor)))
total = sum(observed)
observed
##                 latitiude_factor
## longitude_factor  < 36 >= 36
##          < -120    287  7840
##          >= -120 11406   900
lat_less <- sum(observed[, "< 36"]) / total
lat_greater <- sum(observed[, ">= 36"]) / total
long_less <- sum(observed["< -120", ]) / total
long_greater <- sum(observed[">= -120", ]) / total
lat <- t(as.matrix(c(lat_less, lat_greater)))
long <- as.matrix(c(long_less, long_greater))
expected <- long %*% lat * total
expected
##          [,1]     [,2]
## [1,] 4650.762 3476.238
## [2,] 7042.238 5263.762
chi_sq <- sum((observed - expected)^2/expected)
p_val <- pchisq(chi_sq, df=df_chisq, lower.tail = FALSE)
p_val
## [1] 0

На уровне значимости 0.05 отвергаем гипотезу о независимости широты и долготы домов.  

Посмотрим на распределение домов по географии.

plot(houses$longitude ~ houses$latitude, pch=16, xlab="Latitude", ylab="Longitude", main="Geographical distribution of houses.")

#some_table <- table(data.frame(houses$latitude, houses$longitude))
#hcol <- brewer.pal(9, "YlGnBu")
#heatmap(some_table, col=hcol)

Тест 5

Двустворонний доверительный интервал уровня 0.95 для разности стоимости домов рядом с бухтой и рядом с океаном. Большой объём данных, поэтому используем нормальное распределение.

se <- se_two(house_list$`NEAR BAY`$median_house_value, house_list$`NEAR OCEAN`$median_house_value, type="unknown")
mean_diff <- mymean(house_list$`NEAR BAY`$median_house_value) - mymean(house_list$`NEAR OCEAN`$median_house_value)
delta1 <- mean_diff - abs(qnorm(p = 0.025)) * se
delta2 <- mean_diff + abs(qnorm(p = 0.025)) * se
c(delta1, delta2)
## [1] 10217.27 10256.61

На уровне значимости 0.95 можем сказать, что цена домов рядом с бухтой на $10217 больше, чем домов рядом с океаном.

Тест 6

Правосторонний доверительный интервал уровня 0.95 для разности возраста домов радом с бухтой и рядом с океаном.

se <- se_two(house_list$`NEAR BAY`$housing_median_age, house_list$`NEAR OCEAN`$housing_median_age, type="unknown")
mean_diff <- mymean(house_list$`NEAR BAY`$housing_median_age) - mymean(house_list$`NEAR OCEAN`$housing_median_age)
delta1 <- mean_diff - abs(qnorm(p = 0.05)) * se
c(delta1, Inf)
## [1] 8.275662      Inf

На уровне значимости 0.95 показали, что дома рядом с бухтой на 8 лет старше домов рядом с океаном.

Проверка на выбросы

Построим диаграмму Кливленда стоимостей домов в 5 группах близости к океану, отсортированных по среднему доходу.

by_income <- houses[order(houses$median_income), ]
by_income$ocean_proximity <- factor(by_income$ocean_proximity)
by_income$color[by_income$ocean_proximity == "<1H OCEAN"] <- 1
by_income$color[by_income$ocean_proximity == "INLAND"] <- 2
by_income$color[by_income$ocean_proximity == "NEAR BAY"] <- 3
by_income$color[by_income$ocean_proximity == "ISLAND"] <- 4
by_income$color[by_income$ocean_proximity == "NEAR OCEAN"] <- 5
dotchart(by_income$median_house_value, groups=by_income$ocean_proximity, pch=16, color=by_income$color, xlab="Median house value", ylab="Median income", main="Outliner analysis of house values based on incomes in 5 geographical groups.")

Во всех группах (кроме ISLAND за малым размером) наблюдаем выбросы – точки, которые расположены далеко относительно основной массы. Точки слева от основной массы показывают, что эти дома слишком дёшевы для соответствующего им среднего дохода по сравнению с основной массой, а точки справа – наоборот, дома слишком дороги.  

Посмотрим на boxplot стоимостей домов в зависимости от близости к океану.

boxplot(median_house_value ~ ocean_proximity, data=by_income, pch=16, xlab="Ocean proximity", ylab="Median house value")

В группе INLAND довольно много выбросов над верхним усом. В группе <1H OCEAN тоже есть чуть-чуть таких выбросов.

Проблема множественного тестирования.

Симулируем данные 25 пациентов по 30000 генов в control и condition. Данные по экспрессии генов будем брать из стандартного нормального распределения.

control = list()
samples <- 25
for(i in seq(1:30000)) {
  control[[i]] <- as.vector(rnorm(samples))
}
condition = list()
for(i in seq(1:30000)) {
  condition[[i]] <- as.vector(rnorm(samples))
}

Определим ДЭГи с помощью t-теста.

c = list()
for(i in seq(1:30000)) {
  control_mean <- mean(control[[i]])
  condition_mean <- mean(condition[[i]])
  t_val <- (condition_mean - control_mean) / se_two(x=control[[i]], y=condition[[i]], type = "equal", sx=1, sy=1)
  df_t <- samples + samples - 2
  c[[i]] <- pt(abs(t_val), df=df_t, lower.tail = FALSE) + pt(-abs(t_val), df=df_t)
}
p_vals <- unlist(c)
plot(p_vals, xlab="gene", ylab="p-value", main="P-values without any correction", ylim=c(0, 1.01))
abline(a=0.05, b = 0, col="red")

sum(p_vals < 0.05) / 30000
## [1] 0.0432

Как мы видим, 4% генов неверно определены как ДЭГи при пороге p-value = 0.05.  

Проведём FWER-коррекцию Holm.

p_holm<- p.adjust(p_vals, method="holm")
plot(p_holm, xlab="gene", ylab="p-value", main="P-values with Holm correction", ylim=c(0, 1.01))
abline(a=0.05, b=0, col="red")

sum(p_holm < 0.05) / 30000
## [1] 0

Ни один ген не определён как ДЭг при заданном пороге.  

Проведём FDR-коррекцию Benjamini-Hochberg.

p_bh <- p.adjust(p_vals, method="BH")
plot(p_bh, xlab="gene", ylab="p-value", main="P-values with Benjamini-Hochberg correction", ylim=c(0, 1.01))
abline(a=0.05, b=0, col="red")

sum(p_bh < 0.05) / 30000
## [1] 0

Ни один ген не определён как ДЭГ.

Линейные модели

Небольшое исследование на максимизацию R^2 и минимизацию RMSE. Будем предсказывать цену домов в зависимости от фич.  

houses$NEAR_BAY <- as.numeric(houses$ocean_proximity == "NEAR BAY")
houses$INLAND <- as.numeric(houses$ocean_proximity == "INLAND")
houses$NEAR_OCEAN <- as.numeric(houses$ocean_proximity == "NEAR OCEAN")
houses$ISLAND <- as.numeric(houses$ocean_proximity == "ISLAND")
houses$H_OCEAN <- as.numeric(houses$ocean_proximity == "<1H OCEAN")
numeric_houses <- select(houses, -c("ocean_proximity", "log10_value"))
#Prepare datasets
indices <- caTools::sample.split(numeric_houses$median_house_value, SplitRatio=0.7)
trainData <- numeric_houses[indices, ]
valiData <- numeric_houses[!indices, ]

Для начала построим наивную модель по всем численным фичам (кроме 1-hot векторов).

model0 <- lm(data=trainData, median_house_value ~ (median_income + longitude + latitude + households +
                                                     total_bedrooms + population + total_rooms + housing_median_age))
summary(model0)
## 
## Call:
## lm(formula = median_house_value ~ (median_income + longitude + 
##     latitude + households + total_bedrooms + population + total_rooms + 
##     housing_median_age), data = trainData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -561763  -43823  -10897   30399  503657 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -3.592e+06  7.487e+04 -47.984  < 2e-16 ***
## median_income       4.008e+04  4.026e+02  99.567  < 2e-16 ***
## longitude          -4.297e+04  8.545e+02 -50.283  < 2e-16 ***
## latitude           -4.307e+04  8.087e+02 -53.254  < 2e-16 ***
## households          6.274e+01  9.147e+00   6.859 7.20e-12 ***
## total_bedrooms      1.085e+02  8.368e+00  12.968  < 2e-16 ***
## population         -4.333e+01  1.363e+00 -31.781  < 2e-16 ***
## total_rooms        -7.570e+00  9.630e-01  -7.861 4.09e-15 ***
## housing_median_age  1.166e+03  5.154e+01  22.619  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 69520 on 14492 degrees of freedom
## Multiple R-squared:  0.6405, Adjusted R-squared:  0.6403 
## F-statistic:  3228 on 8 and 14492 DF,  p-value: < 2.2e-16
predData <- predict(model0, valiData)
sqrt(sum((predData - valiData$median_house_value)^2) / dim(valiData)[1])
## [1] 69829.71

Используя линейную комбинацию всех фичей, мы научили модель как-то предсказывать стоимость домов. Однако RMSE высока в целом.  

Добавим 1-hot векторы.

model01 <- lm(data=trainData, median_house_value ~ (median_income + longitude + latitude + households +
                                                     total_bedrooms + population + total_rooms + housing_median_age +
                                                      NEAR_BAY + INLAND + NEAR_OCEAN + ISLAND + H_OCEAN))
summary(model01)
## 
## Call:
## lm(formula = median_house_value ~ (median_income + longitude + 
##     latitude + households + total_bedrooms + population + total_rooms + 
##     housing_median_age + NEAR_BAY + INLAND + NEAR_OCEAN + ISLAND + 
##     H_OCEAN), data = trainData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -555293  -42890   -9934   28926  494483 
## 
## Coefficients: (1 not defined because of singularities)
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -2.241e+06  1.050e+05 -21.331  < 2e-16 ***
## median_income       3.898e+04  4.040e+02  96.483  < 2e-16 ***
## longitude          -2.657e+04  1.218e+03 -21.816  < 2e-16 ***
## latitude           -2.544e+04  1.202e+03 -21.167  < 2e-16 ***
## households          6.286e+01  9.024e+00   6.967 3.39e-12 ***
## total_bedrooms      9.576e+01  8.282e+00  11.563  < 2e-16 ***
## population         -4.281e+01  1.352e+00 -31.671  < 2e-16 ***
## total_rooms        -5.299e+00  9.602e-01  -5.519 3.47e-08 ***
## housing_median_age  1.091e+03  5.204e+01  20.965  < 2e-16 ***
## NEAR_BAY           -5.306e+03  2.256e+03  -2.352   0.0187 *  
## INLAND             -4.073e+04  2.079e+03 -19.596  < 2e-16 ***
## NEAR_OCEAN          4.105e+03  1.865e+03   2.201   0.0277 *  
## ISLAND              1.692e+05  3.962e+04   4.271 1.96e-05 ***
## H_OCEAN                    NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 68550 on 14488 degrees of freedom
## Multiple R-squared:  0.6505, Adjusted R-squared:  0.6503 
## F-statistic:  2248 on 12 and 14488 DF,  p-value: < 2.2e-16
predData <- predict(model01, valiData)
sqrt(sum((predData - valiData$median_house_value)^2) / dim(valiData)[1])
## [1] 69045.17

Сильно лучше не стало.  

Попробуем обойтись минимумом фич. Посмотрим на корреляции фич исходного датасета.

#diverging_col <- rev(brewer.pal(9, "RdBu"))
cor_map <- cor(numeric_houses)
cols <- colorRampPalette(c("blue", "white", "red"))(20)
heatmap(cor_map, col=cols)

Видим явную корреляцию между стоимостью домов и средним доходом. Видим корреляции между фичами (что плохо). Также, из предыдущих заданий, знаем, что цена зависит от географии. Для построения модели 1 выберем средний доход и долготу (широту не будем брать, так как скоррелирована с долготой).

#Test the first model
model1 <- lm(data = trainData, median_house_value ~ (median_income + longitude))
summary(model1)
## 
## Call:
## lm(formula = median_house_value ~ (median_income + longitude), 
##     data = trainData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -546791  -56356  -16486   36785  436649 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -185453.2    41690.3  -4.448 8.72e-06 ***
## median_income   41955.5      367.5 114.152  < 2e-16 ***
## longitude       -1926.2      348.6  -5.526 3.33e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 84050 on 14498 degrees of freedom
## Multiple R-squared:  0.4744, Adjusted R-squared:  0.4743 
## F-statistic:  6542 on 2 and 14498 DF,  p-value: < 2.2e-16
predData <- predict(model1, valiData)
sqrt(sum((predData - valiData$median_house_value)^2) / dim(valiData)[1])
## [1] 82653.61

Модель довольно плохо объясняет данные (R^2 < 0.5).  

Поищем нетривиальные фичи, которые имеют какую-то значимую корреляцию.

plot(houses$median_house_value ~ I((- houses$latitude - houses$longitude)))

cor(houses$median_house_value, -houses$latitude - houses$longitude)
## [1] 0.4912101

Мы видим корреляцию между стоимостью дома и суммой широты и долготы.Построим на этом линейную модель.

model2 <- lm(data=trainData, median_house_value ~ (median_income + I(latitude + longitude)))
summary(model2)
## 
## Call:
## lm(formula = median_house_value ~ (median_income + I(latitude + 
##     longitude)), data = trainData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -506924  -48466  -15244   32236  501455 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             -4069999.3    66251.6  -61.43   <2e-16 ***
## median_income              36688.5      337.9  108.58   <2e-16 ***
## I(latitude + longitude)   -49264.1      793.0  -62.12   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 74770 on 14498 degrees of freedom
## Multiple R-squared:  0.584,  Adjusted R-squared:  0.5839 
## F-statistic: 1.018e+04 on 2 and 14498 DF,  p-value: < 2.2e-16
predData2 <- predict(model2, valiData)
sqrt(sum((predData2 - valiData$median_house_value)^2) / dim(valiData)[1])
## [1] 73818.83
cor(data.frame(predData2, valiData$median_house_value))
##                             predData2 valiData.median_house_value
## predData2                   1.0000000                   0.7633202
## valiData.median_house_value 0.7633202                   1.0000000

В целом, мы повысили R^2 на 0.1 (на 20%), получив 0.58, что уже неплохо. Однако RMSE снизилась несильно (на 10%). Будем искать ещё фичи.

plot(houses$median_house_value ~ I(houses$total_rooms / houses$total_bedrooms))

cor(houses$median_house_value, houses$total_rooms / houses$total_bedrooms)
## [1] 0.3839203

что-то такое есть в этой зависимости. Добавим эту фичу в модель.

model3 <- lm(data=trainData, median_house_value ~ (median_income + I(latitude + longitude) + I((total_rooms / total_bedrooms))))
summary(model3)
## 
## Call:
## lm(formula = median_house_value ~ (median_income + I(latitude + 
##     longitude) + I((total_rooms/total_bedrooms))), data = trainData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -575620  -46803  -14590   31369  513641 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     -3471584.4    70548.2  -49.21   <2e-16 ***
## median_income                      46381.8      549.8   84.37   <2e-16 ***
## I(latitude + longitude)           -42826.7      832.4  -51.45   <2e-16 ***
## I((total_rooms/total_bedrooms))   -19171.4      866.2  -22.13   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 73540 on 14497 degrees of freedom
## Multiple R-squared:  0.5976, Adjusted R-squared:  0.5975 
## F-statistic:  7177 on 3 and 14497 DF,  p-value: < 2.2e-16
predData3 <- predict(model3, valiData)
sqrt(sum((predData3 - valiData$median_house_value)^2) / dim(valiData)[1])
## [1] 72667.61

Чуть-чуть улучшили модель. Видимо, это какой-то потолок для данного подхода.  

Попробуем 1-hot векторы по близости от океана.

model4 <- lm(data=trainData, median_house_value ~ (median_income + NEAR_BAY + INLAND + NEAR_OCEAN + ISLAND + H_OCEAN))
summary(model4)
## 
## Call:
## lm(formula = median_house_value ~ (median_income + NEAR_BAY + 
##     INLAND + NEAR_OCEAN + ISLAND + H_OCEAN), data = trainData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -508758  -46878  -13005   28892  425048 
## 
## Coefficients: (1 not defined because of singularities)
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    84806.7     1702.4  49.817  < 2e-16 ***
## median_income  37016.5      335.4 110.363  < 2e-16 ***
## NEAR_BAY       19540.9     2057.6   9.497  < 2e-16 ***
## INLAND        -79259.6     1479.5 -53.574  < 2e-16 ***
## NEAR_OCEAN     16434.7     1960.2   8.384  < 2e-16 ***
## ISLAND        205541.6    42913.6   4.790 1.69e-06 ***
## H_OCEAN             NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 74300 on 14495 degrees of freedom
## Multiple R-squared:  0.5892, Adjusted R-squared:  0.5891 
## F-statistic:  4159 on 5 and 14495 DF,  p-value: < 2.2e-16
predData4 <- predict(model4, valiData)
sqrt(sum((predData4 - valiData$median_house_value)^2) / dim(valiData)[1])
## [1] 73590.3

В целом, средний доход и близость к океану неплохо объясняют стоимость дома, вполне себе на уровне модели 3 с искусственными фичами. Но модель не объясняет дисперсию внутри групп близости к океану. Поэтому добавим те самые искуственнные фичи и посмотрим, что получится.

model5 <- lm(data=trainData, median_house_value ~ (median_income + I(latitude + longitude) + I(total_rooms / total_bedrooms) +
                                                     NEAR_BAY + INLAND + NEAR_OCEAN + ISLAND + H_OCEAN))
summary(model5)
## 
## Call:
## lm(formula = median_house_value ~ (median_income + I(latitude + 
##     longitude) + I(total_rooms/total_bedrooms) + NEAR_BAY + INLAND + 
##     NEAR_OCEAN + ISLAND + H_OCEAN), data = trainData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -560088  -44888  -12386   30007  485468 
## 
## Coefficients: (1 not defined because of singularities)
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   -1895995.6   104656.0 -18.116  < 2e-16 ***
## median_income                    44415.0      557.7  79.639  < 2e-16 ***
## I(latitude + longitude)         -24059.0     1240.9 -19.389  < 2e-16 ***
## I(total_rooms/total_bedrooms)   -15916.1      877.1 -18.146  < 2e-16 ***
## NEAR_BAY                         16971.6     2019.8   8.403  < 2e-16 ***
## INLAND                          -37995.7     2173.7 -17.480  < 2e-16 ***
## NEAR_OCEAN                       10014.7     1951.0   5.133 2.89e-07 ***
## ISLAND                          178502.7    41843.6   4.266 2.00e-05 ***
## H_OCEAN                               NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 72430 on 14493 degrees of freedom
## Multiple R-squared:  0.6098, Adjusted R-squared:  0.6096 
## F-statistic:  3235 on 7 and 14493 DF,  p-value: < 2.2e-16
predData5 <- predict(model5, valiData)
sqrt(sum((predData5 - valiData$median_house_value)^2) / dim(valiData)[1])
## [1] 71611.55

Модель стала чуть лучше и приблизилась к наивной модели. Посмотрим на корреляции фич между собой.

numeric_houses$geo <- - numeric_houses$latitude - numeric_houses$longitude
numeric_houses$room_ratio <- numeric_houses$total_rooms / numeric_houses$total_bedrooms
cor_map <- cor(numeric_houses)
heatmap(cor_map, col=cols)

total_rooms / total_bedrooms коррелирует со средним доходом, география и стоимость дома антикоррелирует с 1-hot вектором INLAND. В целом, модель 5 вполне неплохо справляется с предсказанием на фоне остальных моделей. Тем не менее, не зря датасет лежит на kaggle :)