Математическое моделирование
Мелихова И.С.
31 марта 2019 г
Данные: Online_Shopping_for_models.csv
library('GGally') # графики совместного разброса переменных## Warning: package 'GGally' was built under R version 3.5.3## Loading required package: ggplot2library('lmtest') # тесты остатков регрессионных моделей## Warning: package 'lmtest' was built under R version 3.5.3## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary('FNN') # алгоритм kNN## Warning: package 'FNN' was built under R version 3.5.3library('mlbench')## Warning: package 'mlbench' was built under R version 3.5.3library('MASS')## Warning: package 'MASS' was built under R version 3.5.3library('ISLR')## Warning: package 'ISLR' was built under R version 3.5.3library('boot')
setwd("D:/Desktop")
DF <- read.table('Online_Shopping_for_models.csv', header = T, # заголовок в первой строке
dec = ',', # разделитель целой и дробной части
sep = ';') # символы пропущенных значений
df <- na.omit(DF)
dim(df)## [1] 11097 18head(df)str(df)## 'data.frame': 11097 obs. of 18 variables:
## $ Administrative : int 1 1 2 11 10 0 2 8 0 0 ...
## $ Administrative_Duration: num 28.2 6 51 471.6 237.8 ...
## $ Informational : int 0 0 0 4 2 0 0 1 0 0 ...
## $ Informational_Duration : num 0 0 0 236 23 ...
## $ ProductRelated : int 1 15 25 22 82 4 16 227 3 1 ...
## $ ProductRelated_Duration: num 0 762 699 883 1815 ...
## $ BounceRates : num 0 0 0 0.00645 0.00233 ...
## $ ExitRates : num 0.05 0.01429 0.00873 0.0182 0.01039 ...
## $ PageValues : num 0 0 0 19.4 8.5 ...
## $ SpecialDay : num 0 0 0 0 0 0.8 0 0 0 0 ...
## $ Month : Factor w/ 10 levels "Aug","Dec","Feb",..: 10 2 7 10 7 7 3 9 7 6 ...
## $ OperatingSystems : int 1 8 2 3 3 2 2 2 2 1 ...
## $ Browser : int 1 13 10 2 2 2 2 2 2 1 ...
## $ Region : int 1 9 1 3 3 1 1 1 3 1 ...
## $ TrafficType : int 5 20 4 4 2 13 6 13 1 9 ...
## $ VisitorType : Factor w/ 3 levels "New_Visitor",..: 1 2 3 3 3 3 3 3 3 3 ...
## $ Weekend : logi FALSE FALSE FALSE FALSE FALSE FALSE ...
## $ Revenue : logi TRUE FALSE FALSE FALSE FALSE FALSE ...my.seed <- 12
set.seed(12)
set.seed(my.seed)
summary(df)## Administrative Administrative_Duration Informational
## Min. : 0.000 Min. : 0.00 Min. : 0.0000
## 1st Qu.: 0.000 1st Qu.: 0.00 1st Qu.: 0.0000
## Median : 1.000 Median : 8.00 Median : 0.0000
## Mean : 2.321 Mean : 80.69 Mean : 0.5047
## 3rd Qu.: 4.000 3rd Qu.: 92.75 3rd Qu.: 0.0000
## Max. :27.000 Max. :3398.75 Max. :24.0000
##
## Informational_Duration ProductRelated ProductRelated_Duration
## Min. : 0.00 Min. : 0.0 Min. : 0.0
## 1st Qu.: 0.00 1st Qu.: 7.0 1st Qu.: 186.5
## Median : 0.00 Median : 18.0 Median : 598.6
## Mean : 34.51 Mean : 31.9 Mean : 1199.3
## 3rd Qu.: 0.00 3rd Qu.: 38.0 3rd Qu.: 1474.5
## Max. :2549.38 Max. :705.0 Max. :63973.5
##
## BounceRates ExitRates PageValues SpecialDay
## Min. :0.000000 Min. :0.00000 Min. : 0.000 Min. :0.00000
## 1st Qu.:0.000000 1st Qu.:0.01423 1st Qu.: 0.000 1st Qu.:0.00000
## Median :0.003139 Median :0.02513 Median : 0.000 Median :0.00000
## Mean :0.022350 Mean :0.04313 Mean : 5.917 Mean :0.06232
## 3rd Qu.:0.016949 3rd Qu.:0.05000 3rd Qu.: 0.000 3rd Qu.:0.00000
## Max. :0.200000 Max. :0.20000 Max. :361.764 Max. :1.00000
##
## Month OperatingSystems Browser Region
## May :3036 Min. :1.000 Min. : 1.000 Min. :1.00
## Nov :2710 1st Qu.:2.000 1st Qu.: 2.000 1st Qu.:1.00
## Mar :1703 Median :2.000 Median : 2.000 Median :3.00
## Dec :1537 Mean :2.123 Mean : 2.352 Mean :3.15
## Oct : 492 3rd Qu.:3.000 3rd Qu.: 2.000 3rd Qu.:4.00
## Sep : 402 Max. :8.000 Max. :13.000 Max. :9.00
## (Other):1217
## TrafficType VisitorType Weekend Revenue
## Min. : 1.000 New_Visitor :1537 Mode :logical Mode :logical
## 1st Qu.: 2.000 Other : 76 FALSE:8512 FALSE:9367
## Median : 2.000 Returning_Visitor:9484 TRUE :2585 TRUE :1730
## Mean : 4.099
## 3rd Qu.: 4.000
## Max. :20.000
## Выдвинем предположение о том, что выручка интернет-магазина будет зависить от показателя отказов, руководства, браузера и соответствующей продукции.
model1 <- lm(Revenue ~ BounceRates + Administrative + Browser + ProductRelated, data = df)
summary(model1)##
## Call:
## lm(formula = Revenue ~ BounceRates + Administrative + Browser +
## ProductRelated, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.85767 -0.17259 -0.13893 -0.08463 1.03913
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.140e-01 6.836e-03 16.671 < 2e-16 ***
## BounceRates -8.427e-01 7.155e-02 -11.778 < 2e-16 ***
## Administrative 8.056e-03 1.136e-03 7.091 1.41e-12 ***
## Browser 6.457e-03 1.961e-03 3.294 0.000992 ***
## ProductRelated 8.425e-04 8.316e-05 10.131 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3545 on 11092 degrees of freedom
## Multiple R-squared: 0.04551, Adjusted R-squared: 0.04517
## F-statistic: 132.2 on 4 and 11092 DF, p-value: < 2.2e-16Все переменные модели являются значимыми. Построим совместный график разброса переменных.
ggp <- ggpairs(df[, c('Revenue','BounceRates', 'Administrative', 'Browser','ProductRelated')],
mapping = ggplot2::aes(color = Revenue))
print(ggp, progress = F)## Warning in ggmatrix_gtable(x, ...): Please use the 'progress' parameter
## in your ggmatrix-like function call. See ?ggmatrix_progress for a few
## examples. ggmatrix_gtable 'progress' and 'progress_format' will soon be
## deprecated.TRUE## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Перекрёстная проверка по отдельным наблюдениям. Это самый затратный в вычислительном плане метод, но и самый надёжный в плане оценки ошибки вне выборки. Попробуем применить его к линейной модели.
fit.glm <- glm(Revenue ~ BounceRates + Administrative + Browser + ProductRelated, data = df)
# считаем LOOCV-ошибку
cv.err <- cv.glm(df, fit.glm)
# результат: первое число -- по формуле LOOCV-ошибки,
# второе -- с поправкой на смещение
cv.err$delta[1]## [1] 0.1257318Ошибка модели оказалась очень низкой. Оценим точность полиномиальных моделей, поочередно меняя степень каждого из регрессора, чтобы найти наилучшую модель.
Регрессор 1
cv.err.loocv <- rep(0, 5)
names(cv.err.loocv) <- 1:5
# цикл по степеням полиномов
for (i in 1:5){
fit.glm <- glm(Revenue ~ poly(BounceRates,i) + Administrative + Browser + ProductRelated, data = df)
cv.err.loocv[i] <- cv.glm(df, fit.glm)$delta[1]
}
# результат
cv.err.loocv## 1 2 3 4 5
## 0.1257318 0.1247703 0.1245650 0.1245298 0.1245363Регрессор 2
cv.err.loocv <- rep(0, 5)
names(cv.err.loocv) <- 1:5
# цикл по степеням полиномов
for (i in 1:5){
fit.glm <- glm(Revenue ~ BounceRates + poly(Administrative,i) + Browser + ProductRelated, data = df)
cv.err.loocv[i] <- cv.glm(df, fit.glm)$delta[1]
}
# результат
cv.err.loocv## 1 2 3 4 5
## 0.1257318 0.1252479 0.1252630 0.1252620 0.1253488Регрессор 3
cv.err.loocv <- rep(0, 5)
names(cv.err.loocv) <- 1:5
# цикл по степеням полиномов
for (i in 1:5){
fit.glm <- glm(Revenue ~ BounceRates + Administrative + poly(Browser,i) + ProductRelated, data = df)
cv.err.loocv[i] <- cv.glm(df, fit.glm)$delta[1]
}
# результат
cv.err.loocv## 1 2 3 4 5
## 0.1257318 0.1257283 0.1257510 0.1257656 0.1257741Регрессор 4
cv.err.loocv <- rep(0, 5)
names(cv.err.loocv) <- 1:5
# цикл по степеням полиномов
for (i in 1:5){
fit.glm <- glm(Revenue ~ BounceRates + Administrative + Browser + poly(ProductRelated,i), data = df)
cv.err.loocv[i] <- cv.glm(df, fit.glm)$delta[1]
}
# результат
cv.err.loocv## 1 2 3 4 5
## 0.1257318 0.1252452 0.1252405 0.1251493 0.1251646Ошибка модели при изменении степени каждого из регрессора практически не меняется. За наилучшую примем первоначальную линейную модель.
Построим график совместного изменения чувствительности и специфичности с изменением вероятности отсечения от 0 до 1 – ROC-кривую. Для примера возьмём модель LDA.
Факт <- df$Revenue
model.logit <- glm(Revenue ~ BounceRates + Administrative + Browser + ProductRelated, data = df, family = 'binomial')
summary(model.logit)##
## Call:
## glm(formula = Revenue ~ BounceRates + Administrative + Browser +
## ProductRelated, family = "binomial", data = df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1891 -0.6274 -0.5755 -0.2276 3.9437
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.792e+00 5.482e-02 -32.695 < 2e-16 ***
## BounceRates -3.040e+01 2.539e+00 -11.975 < 2e-16 ***
## Administrative 4.858e-02 7.774e-03 6.248 4.15e-10 ***
## Browser 4.128e-02 1.470e-02 2.807 0.005 **
## ProductRelated 4.589e-03 5.376e-04 8.537 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 9605.7 on 11096 degrees of freedom
## Residual deviance: 8930.5 on 11092 degrees of freedom
## AIC: 8940.5
##
## Number of Fisher Scoring iterations: 7p.logit <- predict(model.logit, df, type = 'response')
Прогноз <- factor(ifelse(p.logit > 0.5, 2, 1),
levels = c(1, 2),
labels = c('FALSE', 'TRUE'))
# матрица неточностей
conf.m <- table(Факт, Прогноз)
conf.m## Прогноз
## Факт FALSE TRUE
## FALSE 9330 37
## TRUE 1711 19model.lda <- lda(Revenue ~ BounceRates + Administrative + Browser + ProductRelated, data = df)
model.lda## Call:
## lda(Revenue ~ BounceRates + Administrative + Browser + ProductRelated,
## data = df)
##
## Prior probabilities of groups:
## FALSE TRUE
## 0.844102 0.155898
##
## Group means:
## BounceRates Administrative Browser ProductRelated
## FALSE 0.02551689 2.114765 2.330202 28.78659
## TRUE 0.00520270 3.440462 2.469942 48.73642
##
## Coefficients of linear discriminants:
## LD1
## BounceRates -11.14537812
## Administrative 0.10654310
## Browser 0.08539891
## ProductRelated 0.01114260p.lda <- predict(model.lda, df, type = 'response')
Прогноз <- factor(ifelse(p.lda$posterior[, 'TRUE'] > 0.5,
2, 1),
levels = c(1, 2),
labels = c('FALSE', 'TRUE'))
# матрица неточностей
conf.m <- table(Факт, Прогноз)
conf.m## Прогноз
## Факт FALSE TRUE
## FALSE 9292 75
## TRUE 1690 40x <- NULL # для (1 - SPC)
y <- NULL # для TPR
# заготовка под матрицу неточностей
tbl <- as.data.frame(matrix(rep(0, 4), 2, 2))
rownames(tbl) <- c('fact.FALSE', 'fact.TRUE')
colnames(tbl) <- c('predict.FALSE', 'predict.TRUE')
# вектор вероятностей для перебора
p.vector <- seq(0, 1, length = 501)
# цикл по вероятностям отсечения
for (p in p.vector){
# прогноз
Прогноз <- factor(ifelse(p.lda$posterior[, 'TRUE'] > p,
2, 1),
levels = c(1, 2),
labels = c('FALSE', 'TRUE'))
# фрейм со сравнением факта и прогноза
df.compare <- data.frame(Факт = Факт, Прогноз = Прогноз)
# заполняем матрицу неточностей
tbl[1, 1] <- nrow(df.compare[df.compare$Факт == 'FALSE' & df.compare$Прогноз == 'FALSE', ])
tbl[2, 2] <- nrow(df.compare[df.compare$Факт == 'TRUE' & df.compare$Прогноз == 'TRUE', ])
tbl[1, 2] <- nrow(df.compare[df.compare$Факт == 'FALSE' & df.compare$Прогноз == 'TRUE', ])
tbl[2, 1] <- nrow(df.compare[df.compare$Факт == 'TRUE' & df.compare$Прогноз == 'FALSE', ])
# считаем характеристики
TPR <- tbl[2, 2] / sum(tbl[2, 2] + tbl[2, 1])
y <- c(y, TPR)
SPC <- tbl[1, 1] / sum(tbl[1, 1] + tbl[1, 2])
x <- c(x, 1 - SPC)
}
# строим ROC-кривую
par(mar = c(5, 5, 1, 1))
# кривая
plot(x, y, type = 'l', col = 'blue', lwd = 3,
xlab = '(1 - SPC)', ylab = 'TPR',
xlim = c(0, 1), ylim = c(0, 1))
# прямая случайного классификатора
abline(a = 0, b = 1, lty = 3, lwd = 2)
# точка для вероятности 0.5
points(x[p.vector == 0.5], y[p.vector == 0.5], pch = 16)
text(x[p.vector == 0.5], y[p.vector == 0.5], 'p = 0.5', pos = 4)
# точка для вероятности 0.2
points(x[p.vector == 0.2], y[p.vector == 0.2], pch = 16)
text(x[p.vector == 0.2], y[p.vector == 0.2], 'p = 0.2', pos = 4)
Видно, что изменение границы отсечения с 0.5 до 0.2 увеличивает чувствительность модели почти в три раза, в то же время и специфичность значительно ухудшается.