load("E://RDir/df.rda")
df=samsungData
s=266
set.seed(s)
n <- 4000+sample(nrow(df)-4000,1)
set.seed(s)
ds=df[sample(nrow(df), n), ]
require(tree)
require(caret)
str(ds, list.len=0)
## 'data.frame': 6647 obs. of 561 variables:
## [list output truncated]
table(complete.cases(ds))
##
## TRUE
## 6647
Видно, что выборка имеет 6647 наблюдений по 561 переменной. Пропущенных значений нет. Из постановки задачи известно, что две последние переменные - категориальные.
str(ds[,560:561])
## 'data.frame': 6647 obs. of 2 variables:
## $ id : Factor w/ 21 levels "1","3","5","6",..: 17 2 18 9 4 11 17 20 10 11 ...
## $ state: Factor w/ 6 levels "laying","sitting",..: 4 6 4 2 2 1 4 1 1 4 ...
table(ds$id)
##
## 1 3 5 6 7 8 11 14 15 16 17 19 21 22 23 25 26 27
## 303 312 280 288 268 249 292 290 291 336 344 329 370 287 341 371 360 337
## 28 29 30
## 344 310 345
table(ds$state)
##
## laying sitting standing walk walkdown walkup
## 1257 1160 1231 1128 899 972
with(ds,table(id,state))
## state
## id laying sitting standing walk walkdown walkup
## 1 40 44 48 82 43 46
## 3 54 49 57 52 45 55
## 5 47 42 54 53 44 40
## 6 50 53 46 52 46 41
## 7 47 38 49 50 43 41
## 8 47 37 49 44 33 39
## 11 49 50 45 52 43 53
## 14 46 49 50 55 39 51
## 15 61 54 47 50 36 43
## 16 66 62 72 48 42 46
## 17 63 60 74 60 41 46
## 19 78 61 66 50 37 37
## 21 82 78 84 45 40 41
## 22 58 57 54 44 35 39
## 23 70 63 59 56 51 42
## 25 65 61 64 67 55 59
## 26 68 70 68 55 47 52
## 27 66 62 68 53 37 51
## 28 71 66 72 47 40 48
## 29 66 51 54 50 44 45
## 30 63 53 51 63 58 57
Разделить выборку на равные по длине обучающее dft и контрольное dfc подмножества.
index = createDataPartition(y = ds$state,p = 0.5,list=F);
dft=ds[index,] #обучающая выборка
dfc=ds[-index,] #контрольная выборка
Предварительно удалив из таблицы переменную id, с помощью функции tree(state~.,data=dft, split = “gini”, mincut=10) построить дерево классификации tr1. При возникновении ошибки вида «maximum depth reached» следует немного увеличить значение параметра mincut.
min=4 #значение mincut
tr1=tree(state~., data=dft[-560], split="gini", mincut=min)
plot(tr1)
Изучить качество полученного дерева tr1 и определить долю ошибок предсказания на обучающей и контрольной выборке.
summary(tr1)
##
## Classification tree:
## tree(formula = state ~ ., data = dft[-560], split = "gini", mincut = min)
## Variables actually used in tree construction:
## [1] "X65" "X64" "X447" "X536" "X448" "X365" "X522" "X367" "X548" "X449"
## [11] "X288" "X63" "X509" "X439" "X104" "X41" "X292" "X549" "X15" "X12"
## [21] "X3" "X14" "X368" "X4" "X291" "X510" "X369" "X370" "X6" "X39"
## [31] "X293" "X37" "X13" "X277" "X366" "X84" "X1" "X71" "X5" "X458"
## [41] "X38" "X89" "X502" "X61" "X10" "X557" "X149" "X163" "X238" "X28"
## [51] "X94" "X523" "X11" "X81" "X75" "X272" "X30" "X92" "X148" "X194"
## [61] "X32" "X133" "X33" "X160" "X42" "X2" "X290" "X289" "X70" "X23"
## [71] "X44" "X58" "X29" "X173" "X256"
## Number of terminal nodes: 404
## Residual mean deviance: 1.098 = 3208 / 2921
## Misclassification error rate: 0.2379 = 791 / 3325
Мы видим, по каким переменным было осуществлено разбиение дерева. Дерево состоит из Number of terminal nodes узлов, процент ошибки на обучающей выборке равен Misclassification error rate
get_error <- function(tr, df, type_of_df = "Обучающая", error_only=F)
{
p = predict(tr, df); #получаем матрицу вероятностей
len=length(df$state); #длина выборки
predicts = vector(length=len); #вектор с наиболее вероятными состояниями
for (i in 1:len)
predicts[i]=names(which(p[i,]==max(p[i,]))[1]); #вычисляем вектор
if (error_only)
as.numeric(table(predicts==df$state)[1]/len) #выведем только долю ошибочной классификации
else{
print("Ошибка классификации составляет");
print(c("Выборка: ",type_of_df));
print(as.numeric(table(predicts==df$state)[1]/len)); #выведем долю ошибочной классификации
}
}
get_error(tr1, dft)
## [1] "Ошибка классификации составляет"
## [1] "Выборка: " "Обучающая"
## [1] 0.2378947
get_error(tr1, dfc, "Контрольная")
## [1] "Ошибка классификации составляет"
## [1] "Выборка: " "Контрольная"
## [1] 0.3356412
Убедимся, что мы не переобучили дерево, уменьшив параметр mincut до 4 (Задание 3). Для этого построим график зависимости величины ошибки классификации от значения этого параметра.
l = 50;
trainerror = vector(length=l)
controlerror = vector(length=l)
for (i in 0:l-1){
tr = tree(state~., data=dft[-560], split="gini", mincut=(min+i));
trainerror[i+1] = get_error(tr, dft,,T);
controlerror[i+1] = get_error(tr, dfc,,T);
}
mincuts=min:(min+l-1)
plot(trainerror~mincuts, xlim=c((min+l-1),min), ylim=c(min(trainerror),max(controlerror)), type='l', col="green", xlab = "Значение параметра mincut", ylab = "Величина ошибки классификации")
points(controlerror~mincuts, type='l', col="blue")
for (i in 1:l)
points(c(0,max(controlerror))~c((min+l-i),(min+l-i)), type='l')
legend('topright',c("Train","Control"),lty=c(1,1),col=c('green','blue'),ncol=1,bty ="o", bg="white")
С помощью функции cv.tree(tr1, K = 10, method = “misclass”) методом 10-кратного скользящего контроля получить список trees2, содержащий ошибки классификации (свойство dev) для разных размеров дерева (свойство size).
trees2 = cv.tree(tr1, K = 10, method = "misclass")
Определить «на глаз» размер MinTreeSize1 оптимального поддерева дерева tr1, изучив график зависимости ошибки классификации от размера поддерева.
plot(trees2$dev~trees2$size, type='l')
Выше (Задание 4) мы получили график величины ошибки классификации от величины параметра mincut (чем меньше mincut, тем “больше” дерево). Но новый график показывает, что оптимальный размер дерева находится в диапазоне (20,60), что идет вразрез с полученным ранее графиком.
MinTreeSize1 = 20
Значение, определенное “на глаз” = 20
Реализовать с использованием списка trees2 нахождение минимального размера MinTreeSize2 дерева tr1, при котором достигается наименьшая величина ошибки предсказания dev.
Внимание! Функция cv.tree может выполняться достаточно долгое время, которое зависит, в частности, от размера выборки, числа признаков и кратности скользящего контроля.
mindev <- min(trees2$dev)
MinTreeSize2 = min(trees2$size[which(trees2$dev == mindev)])
MinTreeSize2
## [1] 18
Значение, определенное автоматически = 18
С помощью функции cv.tree, указав оптимальный размер (best = MinTreeSize1) дерева, построить это дерево и сохранить в виде объекта tr2.
tr2 = cv.tree(tr1, best= MinTreeSize1)
plot(tr2)
С помощью функции cv.tree, указав оптимальный размер (best = MinTreeSize2) дерева, построить это дерево и сохранить в виде объекта tr3.
tr3 = cv.tree(tr1, best= MinTreeSize2)
plot(tr3)
С помощью функции prune.tree оценить доли ошибок классификации, допускаемых деревьями tr2 и tr3 на обучающей и контрольной выборках.
pt1t=prune.tree(tr2,newdata = dft)
pt1c=prune.tree(tr2,newdata = dfc)
pt2t=prune.tree(tr3,newdata = dft)
pt2c=prune.tree(tr3,newdata = dfc)
plot(pt1t$dev~pt1t$size, type='l', xlab="Size", ylab="Dev")
points(pt1c$dev~pt1c$size, type='l', col='red')
points(pt2t$dev~pt2t$size, type='l', col='green')
points(pt2c$dev~pt2c$size, type='l', col='blue')
legend('topright',c("MinSizeOnEye on Train","MinSizeOnEye on Control", "MinSizeActual on Train", "MinSizeActual on Control"),
lty=c(1,1,1,1),col=c('black','red','green','blue'),ncol=1,bty ="o", bg="white")
errors = vector(length=4)
errors[1]=get_error(tr2,dft,,T)
errors[2]=get_error(tr2,dfc,,T)
errors[3]=get_error(tr3,dft,,T)
errors[4]=get_error(tr3,dfc,,T)
names(errors)=c("tr2 on train","tr2 on control","tr3 on train","tr3 on control")
errors
## tr2 on train tr2 on control tr3 on train tr3 on control
## 0.5392481 0.5439494 0.5557895 0.5571945
Как видно, классификация при помощи поддеревьев имеет неудовлетворительный уровень ошибки.