Данные: Online_Shopping_for_models.csv

Модели: SVM

Расчистаем основные параметры описательной статистики.

library('GGally')       # графики совместного разброса переменных
library('lmtest') # тесты остатков регрессионных моделей
library('FNN') # алгоритм kNN
library('mlbench')
library('ISLR')
library('e1071')     # SVM
library('ROCR')      # ROC-кривые
library('randomForest')      # случайный лес randomForest()
library('gbm')               # бустинг gbm()
library('tree') 
setwd("D:/Desktop")
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  
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)
train <- sample(n, n * train.percent)

Деревья решений

Построим дерево для категориального отклика Revenue.

Revenue <- df$Revenue
Revenue <- as.factor(Revenue)
tree.shopping <- tree(Revenue ~ ., df, subset = train)
summary(tree.shopping)
## 
## Regression tree:
## tree(formula = Revenue ~ ., data = df, subset = train)
## Variables actually used in tree construction:
## [1] "PageValues"              "Administrative_Duration"
## [3] "Month"                   "Informational_Duration" 
## [5] "VisitorType"             "Administrative"         
## [7] "TrafficType"             "Informational"          
## [9] "ProductRelated_Duration"
## Number of terminal nodes:  14 
## Residual mean deviance:  0.05338 = 15.69 / 294 
## Distribution of residuals:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.87500 -0.03226  0.00000  0.00000  0.00000  0.96770
# визуализация
plot(tree.shopping)
text(tree.shopping, pretty = 0)

yhat <- predict(tree.shopping, newdata = df[-train, ])
df.test <- df[-train, "Revenue"]

Рассмотрим более сложные методы улучшения качества дерева. Бэггинг – частный случай случайного леса с \(m = p\), поэтому и то, и другое можно построить функцией randomForest().

Для начала используем бэггинг, причём возьмём все 17 предикторов на каждом шаге (аргумент mtry).

# бэггинг с 14 предикторами
df.test <- df[-train,]
set.seed(my.seed)
bag.df <- randomForest(Revenue ~ ., data = df, subset = train, 
                           mtry = 17, importance = TRUE)
## Warning in randomForest.default(m, y, ...): The response has five or fewer
## unique values. Are you sure you want to do regression?

Можно изменить число деревьев с помощью аргумента ntree.

# бэггинг с 13 предикторами и 25 деревьями
bag.df <- randomForest(Revenue ~ ., data = df, subset = train,
                           mtry = 16, ntree = 25)
## Warning in randomForest.default(m, y, ...): The response has five or fewer
## unique values. Are you sure you want to do regression?

Теперь попробуем вырастить случайный лес. Берём 6 предикторов на каждом шаге.

# обучаем модель
set.seed(my.seed)
rf.df <- randomForest(Revenue ~ ., data = df, subset = train,
                          mtry = 6, importance = TRUE)
## Warning in randomForest.default(m, y, ...): The response has five or fewer
## unique values. Are you sure you want to do regression?
# важность предикторов
importance(rf.df)  # оценки 
##                             %IncMSE IncNodePurity
## Administrative           7.01243560     2.2514137
## Administrative_Duration 15.20309178     3.8685014
## Informational            4.21484408     0.8031463
## Informational_Duration   2.63838479     1.2950842
## ProductRelated           3.01028245     2.5711900
## ProductRelated_Duration  3.16214100     3.0978554
## BounceRates              0.84008803     1.4787549
## ExitRates                3.63764484     3.0654662
## PageValues              31.15494004     8.3153534
## SpecialDay              -0.18903096     0.1761869
## Month                    1.28034160     3.5383491
## OperatingSystems        -0.02029111     0.5559958
## Browser                 -1.40030795     0.4402597
## Region                  -2.29699316     1.0030740
## TrafficType              0.93838479     1.5915002
## VisitorType              6.88040137     0.9903911
## Weekend                  0.80777991     0.4295586
varImpPlot(rf.df)  # графики

По полученным данным можно сделать вывод о том, что наибольшее влияние в модели оказывают такие показатели, как PageValues и Administrative_Duration.

Бустинг

Построим 5000 регрессионных деревьев с глубиной 4.

library(gbm)
df$Weekend <- as.factor(df$Weekend)
set.seed(my.seed)
boost.df <- gbm(Revenue ~ ., data = df[train, ], distribution = "gaussian",
                    n.trees = 5000, interaction.depth = 4)
# график и таблица относительной важности переменных
summary(boost.df)

Теперь построим графики частной зависимости для двух наиболее важных предикторов: PageValues и Month.

par(mfrow = c(1, 2))
plot(boost.df, i = "PageValues")

plot(boost.df, i = "Month")

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

Модель 1: \(\hat{Revenue} = \hat{\beta}_0 + \hat{\beta}_1 \cdot PageValues+\hat{\beta}_2 \cdot Administrative_Duration\).

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

Модель 2: \(\hat{Revenue} = \hat{\beta}_0 + \hat{\beta}_1 \cdot PageValues+\hat{\beta}_2 \cdot Month\).

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

Оба значения MSE оказалась минимальными, с незначительной разницей значение первой модели оказалось меньше, будем использовать эту модель для дальнейшей работы.

Машины опорных векторов

# таблица с данными, отклик — фактор 
PageValues <- df$PageValues
Administrative_Duration <- df$Administrative_Duration
Revenue <- df$Revenue
Revenue <- as.factor(Revenue)
dat <- data.frame(PageValues, Administrative_Duration, Revenue)
# обучающая выборка 
train <- sample(1:nrow(dat), nrow(dat)/2)
# SVM с радиальным ядром и маленьким cost
svmfit <- svm(Revenue ~ ., data = dat[train, ], kernel = "radial", 
              gamma = 1, cost = 1)
plot(svmfit, dat[train, ])

summary(svmfit)
## 
## Call:
## svm(formula = Revenue ~ ., data = dat[train, ], kernel = "radial", 
##     gamma = 1, cost = 1)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1 
##       gamma:  1 
## 
## Number of Support Vectors:  77
## 
##  ( 40 37 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  FALSE TRUE
# SVM с радиальным ядром и большим cost
svmfit <- svm(Revenue ~ ., data = dat[train, ], kernel = "radial", 
              gamma = 1, cost = 1e5)
plot(svmfit, dat[train, ])

summary(svmfit)
## 
## Call:
## svm(formula = Revenue ~ ., data = dat[train, ], kernel = "radial", 
##     gamma = 1, cost = 1e+05)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1e+05 
##       gamma:  1 
## 
## Number of Support Vectors:  64
## 
##  ( 39 25 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  FALSE TRUE

Перекрестная проверка.

# перекрёстная проверка
set.seed(my.seed)
tune.out <- tune(svm, Revenue ~ ., data = dat[train, ], kernel = "radial", 
                 ranges = list(cost = c(0.1, 1, 10, 5),
                               gamma = c(0.5, 0.1,0.05,1, 2, 3)))
summary(tune.out)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost gamma
##     5   0.1
## 
## - best performance: 0.09387097 
## 
## - Detailed performance results:
##    cost gamma      error dispersion
## 1   0.1  0.50 0.12333333 0.05640363
## 2   1.0  0.50 0.10354839 0.06910047
## 3  10.0  0.50 0.11354839 0.05933648
## 4   5.0  0.50 0.10688172 0.06059829
## 5   0.1  0.10 0.12333333 0.05640363
## 6   1.0  0.10 0.10043011 0.05547155
## 7  10.0  0.10 0.09720430 0.06057920
## 8   5.0  0.10 0.09387097 0.06330210
## 9   0.1  0.05 0.12333333 0.05640363
## 10  1.0  0.05 0.10376344 0.05420429
## 11 10.0  0.05 0.09720430 0.06057920
## 12  5.0  0.05 0.10043011 0.05547155
## 13  0.1  1.00 0.12333333 0.05640363
## 14  1.0  1.00 0.11010753 0.06085848
## 15 10.0  1.00 0.11666667 0.06277059
## 16  5.0  1.00 0.12000000 0.06267218
## 17  0.1  2.00 0.12333333 0.05640363
## 18  1.0  2.00 0.11666667 0.06635231
## 19 10.0  2.00 0.11989247 0.06246285
## 20  5.0  2.00 0.11666667 0.05897174
## 21  0.1  3.00 0.12333333 0.05640363
## 22  1.0  3.00 0.11989247 0.06246285
## 23 10.0  3.00 0.11021505 0.04826149
## 24  5.0  3.00 0.11666667 0.06090079
train <- sample(n, n * train.percent)
dat <- data.frame(PageValues, Administrative_Duration, Revenue)
#матрица неточностей для прогноза по лучшей модели
matrix <- table(true = dat[-train, "Revenue"], 
      pred = predict(tune.out$best.model, newdata = dat[-train, ]))
bestmod <- tune.out$best.model
summary(bestmod)
## 
## Call:
## best.tune(method = svm, train.x = Revenue ~ ., data = dat[train, 
##     ], ranges = list(cost = c(0.1, 1, 10, 5), gamma = c(0.5, 
##     0.1, 0.05, 1, 2, 3)), kernel = "radial")
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  5 
##       gamma:  0.1 
## 
## Number of Support Vectors:  71
## 
##  ( 36 35 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  FALSE TRUE
#MSE
sum(diag(matrix))/sum(matrix) 
## [1] 0.8863636

Построим матрицу неточностей для прогноза по лучшей модели на прогнозных данных и рассчитаем MSE.

setwd("D:/Desktop")
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)
Revenue <- df1$Revenue
Revenue <- Revenue[1:30]
PageValues <- df1$PageValues
PageValues <- PageValues[1:30]
Administrative_Duration <- df1$Administrative_Duration
Administrative_Duration <- Administrative_Duration[1:30]
n <- nrow(df1)
# доля обучающей выборки
train.percent <- 0.5
# выбрать наблюдения в обучающую выборку
set.seed(my.seed)
train <- sample(n, n * train.percent)
dat1 <- data.frame(PageValues, Administrative_Duration, Revenue)
#матрица неточностей для прогноза по лучшей модели
matrix <- table(true = dat1[-train, "Revenue"], 
      pred = predict(tune.out$best.model, newdata = dat1[-train, ]))
bestmod <- tune.out$best.model
summary(bestmod)
## 
## Call:
## best.tune(method = svm, train.x = Revenue ~ ., data = dat[train, 
##     ], ranges = list(cost = c(0.1, 1, 10, 5), gamma = c(0.5, 
##     0.1, 0.05, 1, 2, 3)), kernel = "radial")
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  5 
##       gamma:  0.1 
## 
## Number of Support Vectors:  71
## 
##  ( 36 35 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  FALSE TRUE
#MSE
sum(diag(matrix))/sum(matrix) 
## [1] 0.4210526

MSE по лучшей модели на прогнозных данных составило- 0.4210526.

Теперь сравним MSE исходной модели дерева на прогнозных данных с MSE итоговой SVM на прогнозных данных.

df1.test <- df1[-train,]
Revenue <- df$Revenue
Revenue <- as.factor(Revenue)
tree.shopping <- tree(Revenue ~ ., df1, subset = train)
prune.shopping <- prune.misclass(tree.shopping, best = 2)
tree.pred <- predict(prune.shopping, df1.test, type = "class")
# матрица неточностей
Revenue.test <- Revenue[-train]
matrix <- table(tree.pred, Revenue.test)
bestmod <- tune.out$best.model
#MSE
sum(diag(matrix))/sum(matrix)
## [1] 0.6201299

MSE исходной модели дерева на прогнозных данных составило- 0.62, что больше, чем MSE итоговой SVM, а значит лучшая модель выбрана верно.