Практическая работа №6

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

Глушков Егор Александрович, гр. 20.М04-мм

---

• Построить дискриминантную функцию.

• Оценить граничное значение.

• Построить матрицу классификации

• Вычислить точность классификации.

---

Построение модели с помощью встроенной функции lda

library(readxl)
addicts <- read_excel("addicts.xlsx")
# View(addicts)

Исследуем переменные на наличие пропусков. Выделим нужные переменные. В качестве категориальной переменной возьмем curwor – занятость.

X <- data.frame(na.omit(addicts[ , c("sstati", "asi1_med", "asi2_emp", "asid3_dyr", "asi7_psy", "curwor")]))

summary(X)

     sstati         asi1_med         asi2_emp        asid3_dyr         asi7_psy        curwor      
 Min.   :23.00   Min.   :0.0000   Min.   :0.0000   Min.   : 0.500   Min.   :0.00   Min.   :0.0000  
 1st Qu.:43.00   1st Qu.:0.0000   1st Qu.:0.6900   1st Qu.: 2.000   1st Qu.:0.12   1st Qu.:0.0000  
 Median :48.00   Median :0.1900   Median :0.9100   Median : 3.000   Median :0.24   Median :0.0000  
 Mean   :48.52   Mean   :0.2447   Mean   :0.7834   Mean   : 3.605   Mean   :0.26   Mean   :0.2708  
 3rd Qu.:54.00   3rd Qu.:0.4200   3rd Qu.:1.0000   3rd Qu.: 4.000   3rd Qu.:0.36   3rd Qu.:1.0000  
 Max.   :72.00   Max.   :1.0000   Max.   :1.0000   Max.   :24.000   Max.   :0.87   Max.   :1.0000  

summary(as.factor(X$curwor))

  0   1 
202  75 

Применим функцию Линейного Дискриминантного анализа lda с curwor в качестве группирующей переменной

library(MASS)

LDA <- lda(X$curwor ~ ., subset(X, select=-c(curwor))); LDA

Call:
lda(X$curwor ~ ., data = subset(X, select = -c(curwor)))

Prior probabilities of groups:
        0         1 
0.7292419 0.2707581 

Group means:
    sstati  asi1_med  asi2_emp asid3_dyr  asi7_psy
0 48.79703 0.2438119 0.8666337  3.784653 0.2798317
1 47.77333 0.2470667 0.5590667  3.120000 0.2066667

Coefficients of linear discriminants:
                   LD1
sstati    -0.007783931
asi1_med   0.102334038
asi2_emp  -4.073131357
asid3_dyr -0.060416539
asi7_psy  -0.874984879

Веса дискриминантной функции:

alpha.LDA <- LDA$scaling; alpha.LDA

                   LD1
sstati    -0.007783931
asi1_med   0.102334038
asi2_emp  -4.073131357
asid3_dyr -0.060416539
asi7_psy  -0.874984879

Матрица классификации (confusion matrix):

pred <- predict(LDA, X)
tab <- table(real=X$curwor, prediction=pred$class); tab

    prediction
real   0   1
   0 176  26
   1  45  30

Точность классификации (accuracy):

acc <- (tab[1] + tab[4]) / sum(tab); acc

[1] 0.7436823

Заметим, что точность вообще говоря низка: baseline-классификатор с одним ответом о принадлежности к бОльшему классу (тут это 0) для абсолютно всех измерений дал бы 0.7292419.

Соотнесём значения дискриминантной функции, предсказанный класс (т.е. значение curwor – 0 или 1), реальный класс (значение curwor в действительности):

LDA.val <- cbind(pred$x, pred.class=as.numeric(pred$class) - 1, real.class=X$curwor); LDA.val[seq(20), ]

           LD1 pred.class real.class
1   0.33361953          0          1
2   2.64724414          1          1
3  -0.93850975          0          0
4  -0.03465237          0          1
5  -0.29307298          0          0
6   0.86599898          0          0
7  -1.05175811          0          0
8   2.14223046          1          1
9  -0.78318198          0          0
10 -1.21733317          0          0
11 -0.87299796          0          0
12  2.38028413          1          1
13  0.59466102          0          0
14 -0.93944647          0          0
15  0.27721632          0          0
16 -0.85217639          0          0
17 -1.70168728          0          0
18 -0.85045500          0          0
19  2.40385969          1          1
20 -1.05488260          0          0

Построение дискриминантной функции “своими руками”

Ввиду сложности “достать” граничное значение с из встроенной функции lda, самостоятельно реализуем модель линейного дискриминантного анализа для двух классов.
Выборочная ковариационная матрица: \[\Sigma = \frac{1}{n_1 + n_2 - 2} ((n_1 - 1)\Sigma_1 + (n_2 - 1)\Sigma_2)\] где \(\Sigma_1\) и \(\Sigma_2\) – ковариационные матрицы каждой популяции (значения curwor).

Веса дискриминантной функции находятся путём решения СЛАУ \(\Sigma \alpha = \mu_1 - \mu_2\) \[z = \alpha^TX, \quad z_i = \alpha^T \mu_i, \quad i=1,2.\]

W1 <- subset(X, curwor == 0, select=-c(curwor))
W2 <- subset(X, curwor == 1, select=-c(curwor))

n1 <- nrow(W1)
n2 <- nrow(W2)
n <- n1 + n2

mu1 <- colMeans(W1)
mu2 <- colMeans(W2)

S1 <- cov(W1)
S2 <- cov(W2)

# S1 <- 0
# for (i in seq(n1)) {
#   x <- as.matrix(W1[i, ] - mu1)
#   S1 <- S1 + t(x) %*% x
# }
# S1 <- S1 / (n1 - 1)

S <- ((n1 - 1) * S1 + (n2 - 1) * S2) / (n - 2)
# S %*% alpha = mu1 - mu2
alpha <- solve(S, mu1 - mu2); alpha

     sstati    asi1_med    asi2_emp   asid3_dyr    asi7_psy 
 0.01062691 -0.13971017  5.56078793  0.08248287  1.19456136 

Граничное значение \[c = \frac{z_1 + z_2}{2} + \ln {\frac{q_2}{q_1}}\] где \(q_1\) и \(q_2\) – приорные вероятности: \(q_i = \frac{n_i}{n}, \quad i = 1, 2\)

z1 <- t(alpha) %*% mu1
z2 <- t(alpha) %*% mu2

sigma2 <- t(alpha) %*% S %*% alpha

q1 <- n1 / n
q2 <- n2 / n

c <- as.numeric((z1 + z2) / 2 + log(q2 / q1)); c

[1] 4.027395

Проверим модель на исходных данных. Матрица классификации (confusion matrix):

X.red <- subset(X, select=-c(curwor))
mypred <- ifelse(as.matrix(X.red) %*% alpha >= c, 0, 1)
tb <- table(pred=mypred, real=X$curwor); tb

    real
pred   0   1
   0 176  45
   1  26  30

Точность классификации (accuracy)

acc <- (tb[1] + tb[4]) / sum(tb); acc

[1] 0.7436823

Заметим, что “самодельный” LDA сделал прогноз, идентичный прогнозу встроенной функции LDA, несмотря на разные веса alpha.

Проверим идентичность меток классов:

sum(mypred - (as.numeric(pred$class) - 1))

[1] 0

Сумма поэлементного вычитания меток равна 0, то есть метки идентичны.

Проиллюстрируем это, выведя метки классов “самописного” lda, метки встроенного lda, а также настоящие метки [первые два столбца в точности совпадают].

cbind(handmade.pred=as.numeric(mypred), LDA.pred=as.numeric(pred$class) - 1, real.labels=X$curwor)[seq(30),]

      handmade.pred LDA.pred real.labels
 [1,]             0        0           1
 [2,]             1        1           1
 [3,]             0        0           0
 [4,]             0        0           1
 [5,]             0        0           0
 [6,]             0        0           0
 [7,]             0        0           0
 [8,]             1        1           1
 [9,]             0        0           0
[10,]             0        0           0
[11,]             0        0           0
[12,]             1        1           1
[13,]             0        0           0
[14,]             0        0           0
[15,]             0        0           0
[16,]             0        0           0
[17,]             0        0           0
[18,]             0        0           0
[19,]             1        1           1
[20,]             0        0           0
[21,]             0        0           0
[22,]             0        0           0
[23,]             0        0           0
[24,]             0        0           1
[25,]             0        0           0
[26,]             0        0           0
[27,]             1        1           0
[28,]             0        0           0
[29,]             1        1           0
[30,]             0        0           0