Xét dữ liệu (thu từ mẫu) gổm một biến phân nhóm và n biến độc lập (`X_1,X_2, . . . , X_n`). Biến phân nhóm có hai giá trị nên ta chia dữ liệu làm hai nhóm A và B, có vectơ trung bình là `bar mb(x)_A` và `bar mb(x)_B`. Ta giả sử rằng ma trận hiệp phương sai `mb(Sigma)` của hai tổng thể chứa hai nhóm bằng nhau. Trong trường hợp tổng quát, số phân tử `p_A` và `p_B` của hai nhóm khác nhau nên ta có ma trận hiệp phương sai gộp `mb(S)_p` của mẫu.

Như đã trình bày ở trên, với hai nhóm, ta chỉ cần 1 chỉ tiêu phân nhóm. Ta ký hiệu chỉ tiêu này là `Z`. Khi ma trân hiệp phương sai của hai nhóm bằng nhau, ta xem `Z` như tổ hợp tuyến tính của các biến độc lập (linear discriminant analysis hay LDA). Vậy:

`Z=a_1X_1+a_2X_2+...+a_iX_i+...+a_nX_n`(1)

Hay dưới dạng vectơ :

`Z=mb(a)^T mb(X)`(2)

`a_i` chính là các hệ số mà ta phải tìm để `Z` thỏa yêu cầu về phân nhóm.

Xét phần tử u thuộc nhóm A, giá trị `z_u` của phần tử này theo trục `Z` được xác định từ các công thức (1) và (2), nghĩa là:

`z_u=a_1X_(1u)+a_2X_(2u)+...+a_iX_(iu)+...+a_nX_(nu) = mb(a)^T mb(x)_u`(3)

Vậy với `p_A` phần tử của nhóm A, ta có trung bình `bar z_A` và phương sai `s_A^2` của `p_A` giá trị `z_u` ấy.

Tương tự, với phần tử `v` thuộc nhóm B, giá trị `z_v` của phần tử này theo trục `Z` cũng được xác định từ các công thức (1) và (2), nghĩa là:

`z_v=a_1X_(1v)+a_2X_(2v)+...+a_iX_(iv)+...+a_nX_(nv) = mb(a)^T mb(x)_v`(4)

Vậy với `p_B` phần tử của nhóm B, ta có trung bình `bar z_B` và phương sai `s_B^2` của `p_B` giá trị `z_v` ấy.

Ngoài ra, do số phân tử `p_A` và `p_B` của hai nhóm khác nhau nên ta có phương sai gộp `s_p^2` của mẫu.

Khoảng cách tương đối suy rộng `D` (khoảng cách Mahalanobis) giữa hai nhóm được xác định như sau:

`D^2=(bar z_A-bar z_B)^2/s_p^2 = [mb(a)^T (bar mb(x)_A - bar mb(x)_B)]^2/(mb(a)^T mb(S)_p mb(a))`

(5)

Để đạt yêu cầu về phân nhóm, `D` phải có giá trị cực đại. Điều này xảy ra khi:

`mb(a)=c\ mb(S)_p^(-1) (bar mb(x)_A - bar mb(x)_B)`(6)

Trong đó, `c` là một số bất kỳ. Như vậy ta có thể tìm được nhiều vectơ `mb(a)` đáp ứng công thức (6), nhưng tất cả các vectơ này đều song song nhau, nghĩa là phương của trục `Z` là duy nhất.

Khi dữ liệu gồm `K` nhóm (`K>2`), ta sẽ suy rộng sự khảo sát ở trên, nghĩa là ta cần tìm vectơ `mb(a)` sao cho khoảng cách giữa `bar z_1, bar z_2, . . . , bar z_k, . . . , bar z_K` là cực đại.

Trước hết, ta nhận xét rằng công thức (5) có thể viết lại dưới dạng:

`D^2=(bar z_A-bar z_B)^2/s_p^2 = [mb(a)^T (bar mb(x)_A - bar mb(x)_B)]^2/(mb(a)^T mb(S)_p mb(a))= ( mb(a)^T (bar mb(x)_A - bar mb(x)_B) (bar mb(x)_A - bar mb(x)_B)^T mb(a) ) / (mb(a)^T mb(S)_p mb(a))`

(7)

Khi ta mở rộng thành `K` nhóm thì ta thay `(bar mb(x)_A - bar mb(x)_B)(bar mb(x)_A - bar mb(x)_B)^T` bằng ma trận `mb(H)` và `mb(S)_p` bằng ma trận `mb(E)`. Hai ma trận này là sự mở rộng của tổng phương sai nhóm `SS_A` và tổng phương sai sai số `SS_E` (xem thêm phần phân tích phương sai nhiều đáp ứng). Khi ấy, nếu ta gọi `lambda` là tỷ số giữa khoảng cách giữa các nhóm và khoảng cách giữa các phần tử trong một nhóm thì ta có:

`lambda=(mb(a)^T mb(Ha)) / (mb(a)^T mb(Ea))`

(8)

Người ta chứng minh được rằng để đạt yêu cầu phân nhóm, `lambda` là các giá trị riêng của ma trận `mb(E)^(-1) mb(H)`. Ma trận này có `s` giá trị được xếp theo thứ tự giảm dần là `lambda_1,lambda_2, . . . , lambda_s`. Với mỗi giá trị ta có một vectơ `mb(a)`, đó chính là vectơ riêng của ma trận `mb(E)^(-1) mb(H)`.

Trong trường hợp phương sai không đồng nhất, ma trận hiệp phương sai của các nhóm không bằng nhau, mô hình tuyến tính ở công thức (1) không đưa đến kết quả hợp lý, ta phải chuyển sang mô hình phi tuyến. Cho đến nay, mô hình bậc hai (quadratic discriminant analysis hay QDA) được sử dụng phổ biến hơn cả.

Trước hết ta chuyển dữ liệu vào R với tên nbh, sau đó dùng lệnh str để xem qua cấu trúc bảng dữ liệu này. Kết quả thu được là:

> str(nbh)
'data.frame':	90 obs. of  7 variables:
 $ Nhom: Factor w/ 3 levels "CP","HP","NP": 2 2 2 2 2 2 2 2 2 2 ...
 $ WID : num  13.5 15.5 14.5 15.5 14.5 14 15 15 15.5 15.5 ...
 $ CIR : num  57.1 58.4 55.9 58.4 58.4 ...
 $ FBE : num  19.5 21 19 20 20 21 19.5 21 20.5 20.5 ...
 $ EYT : num  12.5 12 10 13.5 13 12 13.5 13 13.5 13 ...
 $ EAT : num  14 16 13 15 15.5 14 15.5 14 14.5 15 ...
 $ JAW : num  11 12 12 12 12 13 13 13 12.5 13 ...

Bảng kết quả trên cho thấy :

• bảng này gồm 90 phần tử,
• biến phụ thuộc là "Nhom" có 3 giá trị tương ứng với 3 nhóm phần tử là CP (college player: vận động viên đại học cộng đồng), HP (highschool player: vận động viên học sinh trung học), NP (non player: không phải là vận động viên),
• 6 biến độc lập tương ứng với 6 kích thước của đầu: WID (chiều ngang đầu), CIR (chu vi đầu), FBE (khoảng cách trước - sau đo tại mắt), EYT (khoảng cách từ mắt đến đỉnh đầu), EAT (khoảng cách từ tai đến đỉnh đầu), JAW (chiều ngang của hàm).

Để có một cái nhìn trực quan hơn về toàn bộ dữ liệu, ta vẽ biểu đồ với lệnh sau:

 plot(nbh, col = nbh$Nhom, pch = 20, cex = 1.1)

Và thu được Hình 3.

Hình 3 Tương quan giữa các biến

Để phân tích sự khác biệt của dữ liệu trên, ta sử dụng lệnh lda (linear discriminant analysis) của phụ kiện MASS. Hai đối số quan trọng của lệnh này là tên bảng dữ liệu và quan hệ giữa biến phụ thuộc và các biến độc lập. Vậy để xử lý dữ liệu trên ta dùng đoạn lệnh sau:

 library(MASS)
 kq <- lda(Nhom ~ ., data = nbh)

Kết quả là một danh sách với 10 thành phần (dùng lệnh names):

> names(kq)
 [1] "prior"   "counts"  "means"   "scaling" "lev"     "svd"     "N"      
 [8] "call"    "terms"   "xlevels"

Để xem xét kết quả, ta dùng lệnh kq và thu được kết quả sau:

> kq
Call:
lda(Nhom ~ ., data = nbh)

Prior probabilities of groups:
       CP        HP        NP 
0.3333333 0.3333333 0.3333333 

Group means:
     WID      CIR      FBE      EYT      EAT      JAW
CP 15.42 57.37967 19.80333 10.08000 13.45333 11.94333
HP 15.20 58.93700 20.10833 13.08333 14.73333 12.26667
NP 15.58 57.77000 19.81000 10.94667 13.69667 11.80333

Coefficients of linear discriminants:
             LD1           LD2
WID -0.948423100  1.4067750094
CIR  0.003639865 -0.0005126312
FBE  0.006439599 -0.0286176430
EYT  0.647483088  0.5402700415
EAT  0.504360916 -0.3839132257
JAW  0.828535064 -1.5288556226

Proportion of trace:
  LD1   LD2 
0.943 0.057

Kết quả trên cho thấy ta có hai chỉ tiêu là LD1 và LD2, mối tương quan giữa hai chỉ tiêu này với các biến độc lập được cho trong phần "Coefficient of linear discriminant". Dựa vào các hệ số trong phần này, ta có thể phán đoán được mối tương quan giữa chỉ tiêu và biến độc lập. Thí dụ chiều ngang đầu (WID) có tương quan chặt chẽ với LD1 và LD2.

Để có cái nhìn trực quan hơn về kết quả phân nhóm bằng các chỉ tiêu LD1 và LD2, ta vẽ biểu đồ xy của các phần tử trên hệ trục LD1 và LD2 này với sự hỗ trợ của phụ kiện ggplot2 và đoạn lệnh sau:

 library(ggplot2)
 kqP <- cbind(scale(as.matrix(nbh[,-1]), scale = FALSE) %*% kq$scaling,
               nbh[,1,drop = FALSE])
 ggplot(data = kqP, aes(x = LD1, y = LD2, col = Nhom, shape = Nhom)) + 
   geom_point() + geom_text(aes(label = Nhom), size = 3)

Trong ba lệnh trên, lệnh thứ hai để tạo một bảng dữ liệu kqP cung cấp cho lệnh thứ ba ggplot. Ta thu được Hình 4.

Hình 4 Kết quả phân nhóm bằng LD1 và LD2

Qua Hình 4, ta thấy LD1 tách riêng nhóm HP, sau đó LD2 tách riêng hai nhóm CP và NP.