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 мало, а значит, распределение не является нормальным. По графику видим, что есть небольшое утяжеление правого хвоста, что, по-видимому, обусловлено распространённостью элитной недвижимости.
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.
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.
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.
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)
Двустворонний доверительный интервал уровня 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 больше, чем домов рядом с океаном.
Правосторонний доверительный интервал уровня 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 :)