Данные: Online_Shopping_for_models.csv

library('GGally')       # графики совместного разброса переменных
## Warning: package 'GGally' was built under R version 3.5.3
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.5.3
library('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.numeric
library('FNN') # алгоритм kNN
## Warning: package 'FNN' was built under R version 3.5.3
library('mlbench')
## Warning: package 'mlbench' was built under R version 3.5.3
library('MASS')
## Warning: package 'MASS' was built under R version 3.5.3
library('ISLR')
## Warning: package 'ISLR' was built under R version 3.5.3
library('boot')
## Warning: package 'boot' was built under R version 3.5.3
DF <- read.table('Online_Shopping_for_models.csv', header = T,            # заголовок в первой строке
                 dec = ',',             # разделитель целой и дробной части
                 sep = ';')     # символы пропущенных значений
df <- na.omit(DF)
dim(df)
## [1] 616  18
head(df)
str(df)
## 'data.frame':    616 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 :   4.00         Median :0.0000  
##  Mean   : 2.239   Mean   :  75.32         Mean   :0.4627  
##  3rd Qu.: 3.000   3rd Qu.:  85.20         3rd Qu.:0.0000  
##  Max.   :24.000   Max.   :2047.23         Max.   :9.0000  
##                                                           
##  Informational_Duration ProductRelated   ProductRelated_Duration
##  Min.   :   0.0         Min.   :  0.00   Min.   :    0.0        
##  1st Qu.:   0.0         1st Qu.:  7.75   1st Qu.:  166.2        
##  Median :   0.0         Median : 16.00   Median :  531.8        
##  Mean   :  26.4         Mean   : 31.08   Mean   : 1158.1        
##  3rd Qu.:   0.0         3rd Qu.: 36.25   3rd Qu.: 1378.3        
##  Max.   :1146.7         Max.   :318.00   Max.   :13717.4        
##                                                                 
##   BounceRates         ExitRates         PageValues       SpecialDay     
##  Min.   :0.000000   Min.   :0.00000   Min.   :  0.00   Min.   :0.00000  
##  1st Qu.:0.000000   1st Qu.:0.01518   1st Qu.:  0.00   1st Qu.:0.00000  
##  Median :0.003689   Median :0.02671   Median :  0.00   Median :0.00000  
##  Mean   :0.024784   Mean   :0.04488   Mean   :  5.66   Mean   :0.05227  
##  3rd Qu.:0.020033   3rd Qu.:0.05000   3rd Qu.:  0.00   3rd Qu.:0.00000  
##  Max.   :0.200000   Max.   :0.20000   Max.   :261.49   Max.   :1.00000  
##                                                                         
##      Month     OperatingSystems    Browser           Region     
##  May    :171   Min.   :1.000    Min.   : 1.000   Min.   :1.000  
##  Nov    :155   1st Qu.:2.000    1st Qu.: 2.000   1st Qu.:1.000  
##  Dec    : 81   Median :2.000    Median : 2.000   Median :3.000  
##  Mar    : 78   Mean   :2.083    Mean   : 2.404   Mean   :3.094  
##  Oct    : 35   3rd Qu.:2.000    3rd Qu.: 2.000   3rd Qu.:4.000  
##  Sep    : 27   Max.   :8.000    Max.   :13.000   Max.   :9.000  
##  (Other): 69                                                    
##   TrafficType                VisitorType   Weekend         Revenue       
##  Min.   : 1.000   New_Visitor      : 85   Mode :logical   Mode :logical  
##  1st Qu.: 2.000   Other            :  3   FALSE:487       FALSE:526      
##  Median : 3.000   Returning_Visitor:528   TRUE :129       TRUE :90       
##  Mean   : 4.244                                                          
##  3rd Qu.: 4.000                                                          
##  Max.   :20.000                                                          
##

Выдвинем предположение о том, что выручка интернет-магазина будет зависить от показателя отказов, руководства, браузера и соответствующей продукции.

Метод проверочной выборки

Он состоит в том, что мы отбираем одну тестовую выборку и будем считать на ней ошибку модели.

# общее число наблюдений
n <- nrow(df)
# доля обучающей выборки
train.percent <- 0.5
# выбрать наблюдения в обучающую выборку
set.seed(my.seed)
inTrain <- sample(n, n * train.percent)

Построим модели для проверки точности. Вид моделей:

Линейная модель: \(\hat{Revenue} = \hat{\beta}_0 + \hat{\beta}_1 \cdot BounceRates+\hat{\beta}_2 \cdot Administrative+\hat{\beta}_3 \cdot Browser+\hat{\beta}_4 \cdot ProductRelated\).

# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(df)
# подгонка линейной модели на обучающей выборке
fit.lm <- lm(Revenue ~ BounceRates +  Administrative + Browser + ProductRelated, 
               subset = inTrain)
# считаем MSE на тестовой выборке
mean((df$Revenue[-inTrain] - predict(fit.lm,
                              df[-inTrain, ]))^2)
## [1] 0.1111203
# отсоединить таблицу с данными
detach(df)

Построим квадратичные модели и оценим MSE при увеличении степеней объясняющих переменных.

Модели:

1)\(\hat{Revenue} = \hat{\beta}_0 + \hat{\beta}_1 \cdot BounceRates^2+\hat{\beta}_2 \cdot Administrative+\hat{\beta}_3 \cdot Browser+\hat{\beta}_4 \cdot ProductRelated\)

2)\(\hat{Revenue} = \hat{\beta}_0 + \hat{\beta}_1 \cdot BounceRates+\hat{\beta}_2 \cdot Administrative^2+\hat{\beta}_3 \cdot Browser+\hat{\beta}_4 \cdot ProductRelated\)

3)\(\hat{Revenue} = \hat{\beta}_0 + \hat{\beta}_1 \cdot BounceRates+\hat{\beta}_2 \cdot Administrative+\hat{\beta}_3 \cdot Browser^2+\hat{\beta}_4 \cdot ProductRelated\)

4)\(\hat{Revenue} = \hat{\beta}_0 + \hat{\beta}_1 \cdot BounceRates+\hat{\beta}_2 \cdot Administrative+\hat{\beta}_3 \cdot Browser+\hat{\beta}_4 \cdot ProductRelated^2\)

Модель 1:

# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(df)
# подгонка линейной модели на обучающей выборке
fit.lm.1 <- lm(Revenue ~ poly(BounceRates,2) +  Administrative + Browser + ProductRelated, 
               subset = inTrain)
# считаем MSE на тестовой выборке
mean((df$Revenue[-inTrain] - predict(fit.lm.1,
                              df[-inTrain, ]))^2)
## [1] 0.1103187
# отсоединить таблицу с данными
detach(df)

Модель 2:

# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(df)
# подгонка линейной модели на обучающей выборке
fit.lm.2 <- lm(Revenue ~ BounceRates +  poly(Administrative,2) + Browser + ProductRelated, 
               subset = inTrain)
# считаем MSE на тестовой выборке
mean((df$Revenue[-inTrain] - predict(fit.lm.2,
                              df[-inTrain, ]))^2)
## [1] 0.1113847
# отсоединить таблицу с данными
detach(df)
## Warning in Auto$mpg[-inTrain] - predict(fit.lm.2, df[-inTrain, ]): длина
## большего объекта не является произведением длины меньшего объекта

Модель 3:

# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(df)
# подгонка линейной модели на обучающей выборке
fit.lm.3 <- lm(Revenue ~ BounceRates +  Administrative + poly(Browser,2) + ProductRelated, 
               subset = inTrain)
# считаем MSE на тестовой выборке
mean((df$Revenue[-inTrain] - predict(fit.lm.3,
                              df[-inTrain, ]))^2)
## [1] 0.1118458
# отсоединить таблицу с данными
detach(df)

Модель 4:

# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(df)
# подгонка линейной модели на обучающей выборке
fit.lm.4 <- lm(Revenue ~ BounceRates +  Administrative + Browser + poly(ProductRelated,2), 
               subset = inTrain)
# считаем MSE на тестовой выборке
mean((df$Revenue[-inTrain] - predict(fit.lm.4,
                              df[-inTrain, ]))^2)
## [1] 0.111589
# отсоединить таблицу с данными
detach(df)

Ошибка модели при изменении степени каждого из регрессора практически не меняется.За наилучшую примем первоначальную линейную модель.Увеличивать степени дальше нет необходимости, так как уже сейчас MSE показывает довольно низкий результат.

Построим модель LDA. Построим матрицу неточностей, оценим чувствительность, специфичность, верность и с помощью этой модели сделаем прогнозы на прогнозных данных.

DF1 <- read.table('Online_Shopping_for_forecast.csv', header = T,            # заголовок в первой строке
                 dec = ',',             # разделитель целой и дробной части
                 sep = ';')     # символы пропущенных значений
df1 <- na.omit(DF1)
Revenue <- c("TRUE", "FALSE")
df1 <- data.frame(df1, Revenue) 
Revenue <- as.factor(Revenue)
Факт <- df$Revenue
Факт <- Факт[1:616]
model.lda <- lda(Revenue ~ BounceRates +  Administrative + Browser + ProductRelated, data = df1)
model.lda
## Call:
## lda(Revenue ~ BounceRates + Administrative + Browser + ProductRelated, 
##     data = df1)
## 
## Prior probabilities of groups:
## FALSE  TRUE 
##   0.5   0.5 
## 
## Group means:
##       BounceRates Administrative  Browser ProductRelated
## FALSE  0.02196897       2.000000 2.340909       29.24675
## TRUE   0.02247122       2.353896 2.545455       29.09740
## 
## Coefficients of linear discriminants:
##                         LD1
## BounceRates     3.928291434
## Administrative  0.251316016
## Browser         0.409015095
## ProductRelated -0.006182088
p.lda <- predict(model.lda, df1, 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   309  217
##   TRUE     56   34
# чувствительность
conf.m[2, 2] / sum(conf.m[2, ])
## [1] 0.3777778
# специфичность
conf.m[1, 1] / sum(conf.m[1, ])
## [1] 0.5874525
# верность
sum(diag(conf.m)) / sum(conf.m)
## [1] 0.5568182

Отчёт по модели LDA содержит три раздела: априарные вероятности классов (Prior probabilities of groups), групповые средние объясняющих переменных (Group means) и коэффициенты линейной разделяющей границы (Coefficients of linear discriminants).

ROC-кривая для LDA

Построим график совместного изменения чувствительности и специфичности с изменением вероятности отсечения от 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  
## -1.5055  -0.5960  -0.5414  -0.2694   2.9232  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -1.717154   0.248217  -6.918 4.58e-12 ***
## BounceRates    -23.543973   9.098154  -2.588  0.00966 ** 
## Administrative   0.077642   0.033467   2.320  0.02034 *  
## Browser         -0.068198   0.073112  -0.933  0.35094    
## ProductRelated   0.004650   0.002484   1.872  0.06126 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 512.38  on 615  degrees of freedom
## Residual deviance: 472.68  on 611  degrees of freedom
## AIC: 482.68
## 
## Number of Fisher Scoring iterations: 7
p.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   525    1
##   TRUE     89    1
model.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.8538961 0.1461039 
## 
## Group means:
##       BounceRates Administrative  Browser ProductRelated
## FALSE 0.027933370       1.984791 2.431559       28.03422
## TRUE  0.006380339       3.722222 2.244444       48.90000
## 
## Coefficients of linear discriminants:
##                        LD1
## BounceRates    -8.45773794
## Administrative  0.15517205
## Browser        -0.08321771
## ProductRelated  0.01027023
p.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   523    3
##   TRUE     85    5
x <- 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)
}

Рассмотрим изменение чувствительности модели при изменение границы отсечения с 0.5 до 0.2.

# строим 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 увеличивает чувствительность модели почти в шесть раз, в то время как специфичность немного ухудшается.

Рассмотрим изменение чувствительности модели при изменение границы отсечения с 0.5 до 0.08.

# строим 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.08
points(x[p.vector == 0.08], y[p.vector == 0.08], pch = 16)
text(x[p.vector == 0.08], y[p.vector == 0.08], 'p = 0.08', pos = 4)

Видно, что изменение границы отсечения с 0.5 до 0.08 увеличивает чувствительность модели почти в восемнадцать раз, но в то же время специфичность ухудшается очень сильно.