Chúng ta sẽ sử dụng dữ liệu do Lebart, Morineau, và Fénelon (1982) cung cấp về chi tiêu của các gia đình Pháp cho thực phẩm. Các dữ liệu này được trình bày ở Bảng 1 (tập tin chi-tieu-thuc-pham.csv).

Bảng 1 Chi tiêu của các gia đình Pháp cho thực phẩm

Gia_Dinh	Thanh_Vien	Banh_Mi	Rau	Trai_Cay	Thit	Gia_Cam	Sua	Ruou_Vang
LD	2	332	428	354	1437	526	247	427
LD	3	406	563	341	1507	544	324	407
LD	4	534	660	367	1620	628	414	407
LD	5	655	776	423	1848	759	495	486
VP	2	293	559	388	1527	567	239	258
VP	3	386	608	396	1501	558	319	363
VP	4	460	699	484	1856	762	400	416
VP	5	584	995	548	2056	893	518	319
QL	2	372	767	562	1948	927	235	433
QL	3	438	843	689	2345	1148	243	341
QL	4	385	789	621	2366	1149	304	282
QL	5	515	1097	887	2630	1167	561	284

Trong Bảng 1 có 9 biến với ý nghĩa như sau:

• Gia_Dinh : loại nghề nghiệp của người chủ gia đình, LD là lao động sử dụng nhiều sức lực, VP là nhân viên văn phòng, QL là cán bộ quản lý cấp trung và cấp cao,
• Thanh_Vien : số thành viên trong gia đình,
• các biến còn lại là chi tiêu cho các loại thực phẩm bánh mì, rau, trái cây, thịt, gia cầm, sữa và rượu vang.

So với nguyên bản, chúng tôi đã sắp xếp lại các dòng theo thứ tự nghề nghiệp và số thành viên, đồng thời tách biến "Family" trong nguyên bản thành hai biến Gia_Dinh và Thanh_Vien để chúng ta có thể khảo sát riêng tác động của hai biến này đến chi tiêu cho thực phẩm.

Trong bảng dữ liệu này, 7 biến Banh_Mi, Rau, Trai_Cay, Thit, Gia_Cam, Sua, Ruou_Vang là các biến phân tích, hai biến Gia_Dinh và Thanh_Vien là các biến phụ.

Trước hết ta nhập bảng dũ liệu vào R và đặt tên là ct. Sau đó xem xét mối tương quan của các biến trong dữ liệu này (ngoại trừ Gia_Dinh) bằng lệnh cor(ct[, -1]) rồi sau đó làm tròn đến 4 chữ số thập phân thì thu được kết quả sau:

> round(cor(ct[, -1]), 4)
           Thanh_Vien Banh_Mi     Rau Trai_Cay    Thit Gia_Cam    Sua Ruou_Vang
Thanh_Vien     1.0000  0.8791  0.7152   0.4027  0.5267  0.4179 0.9275   -0.0504
Banh_Mi        0.8791  1.0000  0.5931   0.1961  0.3213  0.2444 0.8556    0.3038
Rau            0.7152  0.5931  1.0000   0.8563  0.8811  0.8261 0.6628   -0.3565
Trai_Cay       0.4027  0.1961  0.8563   1.0000  0.9595  0.9263 0.3322   -0.4863
Thit           0.5267  0.3213  0.8811   0.9595  1.0000  0.9818 0.3746   -0.4372
Gia_Cam        0.4179  0.2444  0.8261   0.9263  0.9818  1.0000 0.2306   -0.4011
Sua            0.9275  0.8556  0.6628   0.3322  0.3746  0.2306 1.0000    0.0069
Ruou_Vang     -0.0504  0.3038 -0.3565  -0.4863 -0.4372 -0.4011 0.0069    1.0000

Bảng kết quả phân tích trên cho thấy có một số cặp biến có hệ số tương quan tuyến tính `r` rất lớn (trên 0,95). Trường hợp này được biết dưới tên dạng thức Heywood. Dạng thức Heywood sẽ làm một số kết quả có giá trị không bình thường như ta thấy trong phần tiếp.

Trong phân tích thành tố chính, ta chọn lựa số thành tố chính sau khi phân tích. Tuy nhiên trong phân tích yếu tố, ta phải lựa chọn yếu tố trước khi phân tích. Để giúp đỡ chúng ta trong việc lựa chọn này, phụ kiện psych có hàm fa.parallel. Sử dụng hàm này để đánh giá các biến phân tích bằng dòng lệnh:

 fa.parallel(ct[,-c(1,2)])

Ta thu được biểu đồ trên Hình 1.

Hình 1 Biểu đồ giá trị riêng theo số yếu tố

Dựa vào kết quả của biểu đô trên, dựa vào kết quả của phân tích từ hàm fa.parallel cũng như dựa vào một vài khuyến cáo khác (như giá trị riêng lớn hơn 1), ta có thể chọn số yếu tố là 2. Tuy nhiên do ảnh hưởng của dạng thức Heywood nên ta chọn thêm 1 yếu tố nữa là 3 yếu tố.

Để phân tích yếu tố bảng dữ liệu trên, ta dùng hàm fa của phu kiện psych. So sánh với hàm factanal của phụ kiện stats (được kích hoạt sẵn), hàm fa cho phép ta có nhiều tùy chọn hơn trong xử lý số liệu, kết quả thu được phong phú hơn, ngoài ra ta cũng liên kiết được với lệnh plot của phụ kiện psych.

Khi sử dụng hàm fa ta cần lưu ý đến một số đối số sau:

• nfactors : số yếu tố cần xác định,
• covar : ma trận sử dụng để xác định ma trận hệ số tải `mb(L)`. Nếu đối số này bằng FALSE (giá trị mặc định), fa sử dụng ma trận tương quan `mb(R)` để xác định `mb(L)`. Nếu đối só này bằng TRUE, ma trận hiệp phương sai `mb(S)` được sử dụng.
• fm : phương pháp sử dụng để phân tích yếu tố. Giá trị đối số này có thể là:
  • "minres" (mặc định) : sử dụng phương pháp phần dư cực tiểu (minimum residual),
  • "wls" : sử dụng phương pháp bình phương cực tiểu có trọng số (weighted least squares),
  • "gls" : sử dụng phương pháp bình phương cực tiểu có trọng số mở rộng (generaliszed weighted least squares),
  • "pa" : sử dụng phương pháp trục chính, còn gọi là yếu tố chính (principal axis hay principal factor),
  • "lm" : sử dụng phương pháp khả hợp cực đại (maximum likelihood),
  • "minchi" : sử dụng phương pháp cực tiểu giá trị `chi^2` có trọng số của mẫu.
• rotate : phương pháp quay được sử dụng trong xử lý. Ngoài việc không quay (tương ứng với giá trị "none"), ta có hai nhóm phương pháp quay chính sau:
  • nhóm các phương pháp quay trực giao " : "varimax", "quartimax", "bentlerT", "equamax", "varimin", "geominT", và "bifactor",
  • nhóm các phương pháp quay không trực giao : "promax", "oblimin" (mặc định), "simplimax", "bentlerQ", "biquartimin", "geominQ", và "cluster".

Để đơn giản hóa quá trình phân tích, ta sử dụng dòng lệnh sau:

 kq <- fa(ctp[,-c(1,2)], nfactors = 3)

Ta thấy xuất hiện cảnh báo sau :

> kq <- fa(ctp[,-c(1,2)], nfactors = 3)
Warning message:
In fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate,  :
   A Heywood case was detected.  Examine the loadings carefully.

Chương trình xử lý phát hiện dạng thức Heywood trong dữ liệu của ta như ta đã xem xét ở trên.

Kết quả phân tích được lưu vào biến kq có kiểu dữ liệu là danh sách. Để tìm hiểu một cách khái quát cấu trúc của kq, ta có thể dùng lệnh names và có kết quả sau:

> names(kq)
 [1] "residual"     "dof"          "chi"          "nh"           "rms"          "EPVAL"       
 [7] "crms"         "EBIC"         "ESABIC"       "fit"          "fit.off"      "sd"          
[13] "factors"      "complexity"   "n.obs"        "objective"    "criteria"     "STATISTIC"   
[19] "PVAL"         "Call"         "null.model"   "null.dof"     "null.chisq"   "TLI"         
[25] "RMSEA"        "BIC"          "SABIC"        "r.scores"     "R2"           "valid"       
[31] "score.cor"    "weights"      "rotation"     "communality"  "uniquenesses" "values"      
[37] "e.values"     "loadings"     "model"        "fm"           "rot.mat"      "Phi"         
[43] "Structure"    "scores"       "r"            "np.obs"       "fn"

Ta thấy kq bao gồm 47 hạng mục. Để có thể xem xét các nội dung chính, ta dùng lệnh kq hay summary và thu được kết quả sau:

> kq
Factor Analysis using method =  minres
Call: fa(r = ctp[, -c(1, 2)], nfactors = 3)

 Warning: A Heywood case was detected. 
Standardized loadings (pattern matrix) based upon correlation matrix
            MR1   MR2   MR3   h2     u2 com
Banh_Mi    0.10  0.85  0.34 0.95 0.0481 1.3
Rau        0.70  0.46 -0.09 0.95 0.0530 1.8
Trai_Cay   0.89  0.04 -0.24 0.98 0.0224 1.1
Thit       0.98  0.04  0.00 0.99 0.0058 1.0
Gia_Cam    1.06 -0.12  0.12 1.00 0.0040 1.1
Sua       -0.04  1.00 -0.14 0.97 0.0316 1.0
Ruou_Vang -0.33  0.13  0.58 0.56 0.4419 1.7

                       MR1  MR2  MR3
SS loadings           3.65 2.08 0.66
Proportion Var        0.52 0.30 0.09
Cumulative Var        0.52 0.82 0.91
Proportion Explained  0.57 0.33 0.10
Cumulative Proportion 0.57 0.90 1.00

 With factor correlations of 
      MR1  MR2   MR3
MR1  1.00 0.33 -0.27
MR2  0.33 1.00  0.12
MR3 -0.27 0.12  1.00

Mean item complexity =  1.3
Test of the hypothesis that 3 factors are sufficient.

The degrees of freedom for the null model are  21  and the objective function was  
14.52 with Chi Square of  113.7
The degrees of freedom for the model are 3  and the objective function was  2.23 

The root mean square of the residuals (RMSR) is  0.02 
The df corrected root mean square of the residuals is  0.05 

The harmonic number of observations is  12 with the empirical chi square  0.16  
with prob <  0.98 
The total number of observations was  12  with MLE Chi Square =  13  with prob <  0.0046 

Tucker Lewis Index of factoring reliability =  -0.099
RMSEA index =  0.807  and the 90 % confidence intervals are  0.256 0.836
BIC =  5.54
Fit based upon off diagonal values = 1
Measures of factor score adequacy             
                                               MR1  MR2  MR3
Correlation of scores with factors               1 0.99 0.96
Multiple R square of scores with factors         1 0.98 0.92
Minimum correlation of possible factor scores    1 0.97 0.85

Ta thấy trên bảng kết quả này xuất hiện gần như đầy đủ các thông tin chính về 3 yếu tố được thể hiện qua các điểm sau:

• Các yếu tố được ký hiệu là MR1, MR2 và MR3.
• Trong phần "Standardized loadings . . .", ta thấy yếu tố thứ nhất MR1 tương quan mật thiết với các biến Rau, Trai_Cay, Thit, Gia_Cam, yếu tố thứ hai MR2 tương quan chặt chẽ với Banh_Mi và Sua, yếu tố thứ ba MR3 có tương quan đáng kể với Ruou_Vang.
• Trong phần đánh giá về sự liên hệ của các yếu tố đến sự biến động của giá trị các biến, MR1 liên hệ đến 57%, MR2 đến 33%, và MR3 đến 10% phương sai của dữ liệu.

Ngoài ra ta còn tìm thấy một số thông tin khác liên quan đến các yếu tố cũng việc đánh giá mô hình như phương sai chung `h^2`, phương sai đặc thù `u^2`, giá trị `chi^2`, giá trị `p`.

Biểu đồ cấu trúc (path diagram) giới thiệu cho ta mối quan hệ giữa biến và yếu tố bằng các đường nối giữa chúng. Để thực hiện biểu đồ này ta sử dụng hàm fa.diagram và lưu ý đến các đối số sau:

• cut : giá trị thấp nhất của hệ số tương quan để được thể hiện trên biểu đồ,
• simple : nếu giá trị đối số này là TRUE (mặc định), chỉ những hệ số tương quan lớn nhất được thể hiện trên biểu đồ; nếu giá trị này là FALSE, tất cả các mối tương quan đều được thể hiện trên biểu đồ,
• digits : số chữ số thập phân của hệ số tương quan.

Trên Hình 2 là biểu đồ cấu trúc của thí dụ trên. Hình 2a tương ứng với lệnh:

 fa.diagram(kq, main = "Phân tích yếu tố")

Còn Hình 2b tương ứng với lệnh :

 fa.diagram(kq, simple = FALSE,  main = "Phân tích yếu tố", digits = 2,
           marg = c(0.5,0.5,1.5,0))

Hình 2 Biểu đồ cấu trúc chi tiêu cho thực phẩm ở Pháp

Biểu đồ tương quan giữa biến và yếu tố, gọi tắt là biểu đồ tương quan, có hai trục thể hiện hai yếu tố, các điểm biểu diễn các biến. Để vẽ biểu đồ này ta sử dụng lệnh fa.plot của phụ kiện psych. Nếu ta chọn đối số choose = NULL (giá trị mặc định), kết quả là một ma trận các biểu đồ tương quan gồm tất cả những sự kết hợp của các yếu tố. Thí dụ như lệnh sau:

 fa.plot(kq, title = "Phân tích tương quan yếu tố - biến",
         labels = colnames(ct[,-c(1, 2)]), xlim = c(-0.6, 1.2), ylim = c(-0.5, 1.2))

Trong lệnh trên, đối số labels = colnames(ct[,-c(1, 2)] dùng để đưa tên các biến vào biểu đồ, các lệnh khác để trang trí, làm biểu đồ cân đối. Kết quả thu được là Hình 3.

Hình 3 Biểu đồ tương quan giữa yếu tố và biến

Trong Hình 3, ta có 6 biểu đồ tùy theo sự kết hợp giữa 3 yếu tố với sự bố trí các trục thành trục hoành và trục tung. Vị trí của các biến cho ta biến mối quan hệ giữa biến và yếu tố. Màu của tên biến và điểm cho ta biết mối quan hệ giữa các biến với nhau. Từ biểu đồ này, ta có một số nhận xét sau:

• Các biến Rau, Trai_Cay, Thit, Gia_Cam tương ứng với giá trị cao của yếu tố MR1. Ngược lại Ruou_Vang lại tương ứng với giá trị thấp của yếu tố này.
• Các biến Banh_Mi và Sua tương ứng với giá trị cao của yếu tố MR2

Nếu ta chỉ muốn vẽ một biểu đồ tương quan ứng với hai yếu tố nào đó thì ta đưa chỉ số của hai yếu tố đó vào đối số choose của lệnh fa.plot. Thí dụ đoạn lệnh sau:

 fa.plot(kq, title = "Phân tích tương quan yếu tố - biến",
        choose = c(1, 2), cex = 0.8 , labels = colnames(ct[,-c(1, 2)]),
        xlim = c(-0.4, 1.2), ylim = c(-0.2, 1.2))

dùng để vẽ biểu đồ tương quan với trục hoành tương ứng với yếu tố MR1, trục tung tương ứng với yếu tố MR2. Kết quả được trình bày trên Hình 4.

Hình 4 Biểu đồ tương quan giữa hai yếu tố MR1 & MR2 và các biến

Ta có thể biểu diễn các phần tử và biến trên cùng một biểu đồ bằng lệnh biplot. Khi ấy ta có thể dùng đối số labels để đưa tên hay một tính chất nào đó của phần tử lên biểu đồ này.

Thí dụ ta muốn đưa nghề nghiệp của chủ gia đình (cột thứ 1 của bảng dữ liệu) lên biểu đồ kép thì ta dùng lệnh sau:

 biplot(kq, choose = c(1, 2), labels = ct[, 1], col = c("blue", "red"),
       xlim.s = c(-1.5, 2.2), ylim.s = c(-1.5, 2), 
       xlim.f = c(-0.6, 1.5), ylim.f = c(-0.2, 1.2))

Và ta thu được biểu đồ trên Hình 5

Hình 5 Biểu đồ kép của phân tích yếu tố

Hay nếu ta muốn đưa số thành viên của gia đình (cột thứ 2 của bảng dữ liệu) lên biểu đồ kép thì ta dùng lệnh sau:

 biplot(kq, choose = c(1, 2), labels = ct[, 2], col = c("blue", "red"),
       xlim.s = c(-1.5, 2.2), ylim.s = c(-1.5, 2), 
       xlim.f = c(-0.6, 1.5), ylim.f = c(-0.2, 1.2))

Và ta thu được biểu đồ trên Hình 6

Hình 6 Biểu đồ kép của phân tích yếu tố

Qua Hình 5 và 6, ta có nhận xét rằng nhóm "QL" có xu hướng tập trung về bên phải của biểu đồ, ứng với giá trị cao của MR1. Mặt khác các gia đình có đông thành viên (5) có xu hướng tập trung ở bên trên của biểu đồ ứng với giá trị cao của MR2.