Классификация неоднородных данных

Кластерный анализ

Реальные данные

Данные о гейзере: faithful

faithful - встроенная выборка данных из библиотеки datasets.
Содержит 272 наблюдения об извержениях большого гейзера в Иеллоустонском национальном парке (США), где переменная eruptions показывает длительность извержений (мин), а переменная waiting - время, которое прошло до следующего извержения (мин).
Справку по данным можно найти с помощью команды ?faithful

#изобразим данные на двухмерном графике
plot(faithful,main="Old Faithful Geyser")

#применим алгоритм кластерного анализа "K-средних" (kmeans)
#алгоритм требует явного задания числа кластеров (классов), зададим 2
res=kmeans(faithful,2)
#просмотрим полученное распределение наблюдений по классам
table(res$cluster)

## 
##   1   2 
## 172 100

#отметим на графике принадлежность к классу
#синий цвет - первый класс, зеленый - второй
plot(faithful,col=ifelse(res$cluster==1,"blue","green"))
#добавим легенду
legend("topleft",legend=c("1","2"),fill=c("blue","green"))

График с использованием различных символов (pch)

plot(faithful,pch=ifelse(res$cluster==1,1,2))
#добавим легенду
legend("topleft",legend=c("1","2"),pch=c(1,2))

Модельные данные

Описание данных

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

Моделирование

#для воспроизводимимости результатов
set.seed(929)

#оценим параметры
#для первого класса
mean1=sapply(faithful[res$cluster==1,],mean)
cov1=cov(faithful[res$cluster==1,])
#для второго класса
mean2=sapply(faithful[res$cluster==2,],mean)
cov2=cov(faithful[res$cluster==2,])
#загрузка библиотеки моделирования многомерных нормальных СВ
library(mnormt)
#смоделируем смесь из многом. норм. СВ
n1=sum(res$cluster==1) #число наблюдений из первого класса
n2=sum(res$cluster==2) #число наблюдений из второго класса
# c помощью функции cbind соединим наблюдения из двух классов в одну таблицу
moddata=rbind(rmnorm(n1,mean1,cov1),rmnorm(n2,mean2,cov2))
colnames(moddata)=c("v1","v2")

#визуализируем полученные данные
plot(moddata)

#создадим вектор истинной классификации
cl=rep(c(1,2),c(n1,n2))

Классификация

#применим алгоритм "K-средних" к модельным данным
res2=kmeans(moddata,2)
cl.est=res2$cluster
#перенумерация классов: установить реальное соответствие номеров классов
r=(nerr=sum(cl.est!=cl))/(n1+n2)
if (r>0.5) cl.est=3-cl.est

Ошибки классификации

#безусловная ошибка классификации: число неверно классифицированных наблюдений
(nerr=sum(cl.est!=cl))

## [1] 7

nerr/(n1+n2)

## [1] 0.02574

#условная ош. кл. для первого класса
(nerr1=sum(cl.est[cl==1]!=1))

## [1] 2

nerr1/n1

## [1] 0.01163

#условная ош. кл. для второго класса
(nerr2=sum(cl.est[cl==2]!=2))

## [1] 5

nerr2/n2

## [1] 0.05

#индексы наблюдений, по которым была допущена ошибка
idx=which(cl.est!=cl)
#обозначим неправильно классифицированные наблюдения на графике красным цветом
plot(moddata,col=ifelse(cl==1,"blue","green"))
legend("topleft",legend=c("1","2"),fill=c("blue","green"))
points(moddata[idx,],col="red")

Дискриминантный анализ

library(MASS)

Формирование тренировочной и тестовой подвыборок

n=nrow(moddata)
#количество наблюдений для обучения: 70%
n.train=floor(n*0.7)
#количество наблюдений для тестирования: 30%
n.test=n-n.train
#выбрать случайным образом n.train наблюдений
idx.train=sample(1:n,n.train) #случайный выборк индексов
data.train=moddata[idx.train,]
#выбрать оставшиеся наблюдения в тестовую выборку
idx.test=(1:n)[!(1:n %in% idx.train)]
data.test=moddata[idx.test,]
#истинные классификации
cl.train=cl[idx.train]
cl.test=cl[idx.test]

Обучение

#оценить модель с помощью обучающей выборки data.train
#используется квадратичный дискриминантный анализ: ковар. матр. различны
#указывается тренировачная выборка и вектор истинной классификации
mod=qda(data.train,cl.train) #результат - оцененная модель

Классификация тестовых данных

#оценить модель с помощью обучающей выборки data.train
#используется квадратичный дискриминантный анализ: ковар. матр. различны
#указывается тренировачная выборка и вектор истинной классификации
cl.test_est=predict(mod,data.test)$class

#оценим безусловную ошибку классификации
sum(cl.test_est!=cl.test)/n.test

## [1] 0.0122

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

library(rpart)
options(digits=7)

Описание данных: kyphosis

Данные о наличии искривления позвоночника (kyphosis) у 81 детей после операции на спинном мозге.  
   Переменные:
1. Kyphosis - номинальная переменная (фактор), характеризующая наличие искревления ("absent" - искривления нет, "present" - есть);
2. Age - возрасть детей (в месяцах);  
3. Number - количество искривленных позвонков;
4. Start - номер самого верхнего из позвонков, которые были оперированы.

Данные

#предварительный просмотр
head(kyphosis)

##   Kyphosis Age Number Start
## 1   absent  71      3     5
## 2   absent 158      3    14
## 3  present 128      4     5
## 4   absent   2      5     1
## 5   absent   1      4    15
## 6   absent   1      2    16

#типы переменных
str(kyphosis)

## 'data.frame':    81 obs. of  4 variables:
##  $ Kyphosis: Factor w/ 2 levels "absent","present": 1 1 2 1 1 1 1 1 1 2 ...
##  $ Age     : int  71 158 128 2 1 1 61 37 113 59 ...
##  $ Number  : int  3 3 4 5 4 2 2 3 2 6 ...
##  $ Start   : int  5 14 5 1 15 16 17 16 16 12 ...

Оценивание модели

#оценить модель
fit <- rpart(Kyphosis ~ Age + Number + Start,
    method="class", data=kyphosis)

#построить дерево решений
plot(fit, uniform=TRUE, 
    main="Classification Tree for Kyphosis")
text(fit, use.n=TRUE, all=TRUE, cex=.8)

Пояснения к графику

1. При выполнении промежуточных условий обход по дереву осуществляется влево, иначе - вправо.
   
2. Числа через дробь означают количество случаев без искривлений и с искривлениями соответственно. При наличии занчительного числа искривлений выводится заголовок “present”, иначе - “absent” (при незначительном их числе.