Introduction

Post trước giới thiệu về SEM (Structural Equation Modeling), ứng dụng trong nghiên cứu của SEM, sơ lược về các tools hiện có được sử dụng cho kĩ thuật phân tích này cũng như trình bày việc tái lập lại kết quả của một case được trình bày trong chương 7 cuốn Marketing Research With SPSS với ngôn ngữ R. Một số phân tích EFA (Exploratory Factor Analysis) chủ yếu (như KMO Test, Cronbach Alpha, Ma trận xoay nhân tố) được trình bày tại đây.

Post này sẽ hướng dẫn việc thực hiện CFA/SEM (Confirmatory Factor Analysis/Structural Equation Modeling) với ngôn ngữ R cũng như việc tính toán và trình bày các kết quả chính trong các nghiên cứu sử dụng hai cases từ cuốn Applied Structural Equation Modeling using AMOS: Basic to Advanced Techniques của Joel E. Collier.

Hướng dẫn download sách tại đây còn data thực hành lấy tại đây.

First Order CFA

Sử dụng bộ dữ liệu Customer Delight Data_Student.sav cho mô hình CFA dưới đây:

Trước hết đọc bộ dữ liệu này vào R:

# Clear R environment: 
rm(list = ls())

# Load data: 
haven::read_sav("Customer Delight Data_Student.sav") -> data_f4.7

Kí hiệu ADA (viết tắt của Adaptive Behavior) là nhân tố cấu thành từ các items/questions có tên adapt1, adapt2, adapt3, adapt4 và adapt5 (đây là các biến theo thang đo Likert 5 hoặc Likert 7). Tương tự với hai nhân tố còn lại kí hiệu là DEL (Customer Delight) và WOM (Word of Mouth). Mô hình CFA trên được tuyên bố như sau:

# Identify our SEM Model: 
my_SEM <- "ADA =~ adapt1 + adapt2 + adapt3 + adapt4 + adapt5
           DEL =~ delight1 + delight2 + delight3
           WOM =~ WOM1 +  WOM2 + WOM3"

Để chạy một hình CFA trên chúng ta sử dụng hàm cfa() (hoặc sem()) của gói lavaan:

# Run CFA model: 
library(lavaan)
SEM_model <- cfa(my_SEM, data = data_f4.7)

Xem qua một số kết quả chính của mô hình:

# Show results: 
summary(SEM_model)
## lavaan 0.6-9 ended normally after 41 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        25
##                                                       
##   Number of observations                           500
##                                                       
## Model Test User Model:
##                                                       
##   Test statistic                                61.041
##   Degrees of freedom                                41
##   P-value (Chi-square)                           0.023
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   ADA =~                                              
##     adapt1            1.000                           
##     adapt2            0.976    0.030   32.773    0.000
##     adapt3            0.998    0.035   28.842    0.000
##     adapt4            0.995    0.034   29.028    0.000
##     adapt5            0.929    0.032   28.607    0.000
##   DEL =~                                              
##     delight1          1.000                           
##     delight2          1.221    0.046   26.568    0.000
##     delight3          1.193    0.050   23.673    0.000
##   WOM =~                                              
##     WOM1              1.000                           
##     WOM2              1.075    0.042   25.769    0.000
##     WOM3              1.061    0.046   23.097    0.000
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   ADA ~~                                              
##     DEL               0.446    0.040   11.249    0.000
##     WOM               0.381    0.053    7.154    0.000
##   DEL ~~                                              
##     WOM               0.379    0.044    8.688    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .adapt1            0.208    0.017   12.496    0.000
##    .adapt2            0.145    0.013   11.285    0.000
##    .adapt3            0.264    0.020   13.224    0.000
##    .adapt4            0.256    0.019   13.160    0.000
##    .adapt5            0.236    0.018   13.302    0.000
##    .delight1          0.178    0.015   11.579    0.000
##    .delight2          0.141    0.017    8.133    0.000
##    .delight3          0.284    0.023   12.105    0.000
##    .WOM1              0.321    0.036    8.916    0.000
##    .WOM2              0.375    0.042    8.979    0.000
##    .WOM3              0.681    0.055   12.340    0.000
##     ADA               0.849    0.066   12.798    0.000
##     DEL               0.491    0.042   11.665    0.000
##     WOM               1.228    0.100   12.234    0.000

Một số chỉ số đánh giá mức độ phù hợp của mô hình CFA:

# Show some goodness of fit indices: 
fitMeasures(SEM_model, 
            fit.measures = c("cfi", "tli", "rmsea", "gfi", "chisq", "df", "fmin", "ifi"))
##    cfi    tli  rmsea    gfi  chisq     df   fmin    ifi 
##  0.996  0.994  0.031  0.979 61.041 41.000  0.061  0.996

Nếu bạn sử dụng AMOS-SPSS thì tất cả các kết quả trên có thể được tái lập lại như sau:

Model Diagram cho mô hình CFA được vẽ bằng hàm lavaanPlot() của thư viện lavaanPlot:

library(lavaanPlot) # Install guide: https://www.alexlishinski.com/post/lavaanplot-0-5-1/

# Show model diagram: 
lavaanPlot(model = SEM_model, coefs = TRUE, covs = TRUE, stars = c("regress"))

Các kết quả của mô hình có thể được trình bày chi tiết hơn bằng hàm model_parameters() của gói parameters:

library(parameters) # For using this package: https://rdrr.io/cran/parameters/man/model_parameters.lavaan.html

# Extract, for example, standardized loadings: 
model_parameters(
  SEM_model,
  standardize = TRUE,
  component = c("loading")) -> loadings_cfa1


# Show standardized loadings: 
loadings_cfa1
## # Loading
## 
## Link            | Coefficient |       SE |       95% CI |      z |      p
## -------------------------------------------------------------------------
## ADA =~ adapt1   |        0.90 |     0.01 | [0.88, 0.92] |  86.83 | < .001
## ADA =~ adapt2   |        0.92 | 8.53e-03 | [0.90, 0.94] | 107.88 | < .001
## ADA =~ adapt3   |        0.87 |     0.01 | [0.85, 0.90] |  72.43 | < .001
## ADA =~ adapt4   |        0.88 |     0.01 | [0.85, 0.90] |  73.76 | < .001
## ADA =~ adapt5   |        0.87 |     0.01 | [0.85, 0.89] |  70.78 | < .001
## DEL =~ delight1 |        0.86 |     0.01 | [0.83, 0.89] |  57.19 | < .001
## DEL =~ delight2 |        0.92 |     0.01 | [0.89, 0.94] |  77.90 | < .001
## DEL =~ delight3 |        0.84 |     0.02 | [0.81, 0.87] |  53.22 | < .001
## WOM =~ WOM1     |        0.89 |     0.01 | [0.86, 0.92] |  62.57 | < .001
## WOM =~ WOM2     |        0.89 |     0.01 | [0.86, 0.92] |  62.30 | < .001
## WOM =~ WOM3     |        0.82 |     0.02 | [0.78, 0.85] |  45.63 | < .001

Report CFA Results

Các mô hình CFA/SEM có rất nhiều outputs bất kể là bạn sử dụng AMOS hay R. Do vậy không phải tất cả các outputs đó sẽ được trình bày trong research paper. Về vấn đề này tác giả J. Collier (p.89) viết:

“There are many acceptable formats to report CFA results, but these templates allow the reader to view all the necessary information and at the same time present it in a succinct manner. Before figuring out the format, you need to assess what information needs to be included. A table should include the actual wording of your indicators as they were presented in the survey; factor loadings, t-values, reliabilities, and model fit statistics are advisable when you want a reader to consider your findings. I include the actual indicator wording in the analysis because it gives a clue to strong or weak indicators. It allows the reader to determine the face validity of your indicators and also can provide insight for the reader on how exactly the construct was captured. Lastly, you need to denote what indicators were constrained for identification purposes or, specifically, which indicator was used to set the metric for each construct.”

Cụ thể, với mô hình CFA trên thì báo cáo cuối cùng cho kết quả sẽ được trình bày như ở Table 4.1 dưới đây:

Các giá trị Model Fit Statistics chính là các chỉ số đánh giá mức độ phù hợp của mô hình và các giá trị t-value thì nằm ở cột z-value của mục Latent Variables khi thực hiện lệnh summary(SEM_model). Các giá trị Loading cho, ví dụ, nhân tố ADA làm tròn ở 2 chữ số lần lượt là 0.89, 0.92, 0.87, 0.87 và 0.86 chúng ta có thể tìm thấy ở loadings_cfa1.

Riêng giá trị C.R = 0.94 cho nhân tố này là Composite Reliability. Theo Yang và Green (2011) thì Cronbach Alpha là một tiêu chí phản ánh không chính xác sự phù hợp (Reliability) của thang đo. Vì lí do này Hair et al. (2009) cho rằng nên thay thế Cronbach Alpha bằng một chỉ số khác là Composite Reliability (viết tắt C.R). Chỉ số này được tính toán dựa trên hệ số tải chuẩn hóa (Standardized Loadings). Ví dụ, với nhân tố ADA thì C.R được tính như sau:

C.R = 0.948 nếu làm tròn 2 chữ số sau dấu phẩy thì sẽ là 0.95 - sai lệch một chút so với 0.94 được trình bày trong Table 4.1. Rất tiếc là AMOS không hỗ trợ tính Composite Reliability nên chúng ta buộc phải tính tay chỉ số này như sau:

# Extract Standardized Loadings for all factors: 
loadings_cfa1$Coefficient -> all_loadings

# Loadings for ADA
all_loadings[1:5] -> ada_loadings

# Calcualte Composite Reliability for ADA: 

ts <- sum(ada_loadings)*sum(ada_loadings)

ms <- (ts + sum(1 - ada_loadings^2))

compositeReliability_ADA <- ts / ms

# Show result: 
compositeReliability_ADA
## [1] 0.9486553

Kết quả không làm tròn là 0.9486553. Chúng ta thực hiện việc tính toán Composite Reliability tương tự cho hai nhân tố còn lại. Nếu quen thuộc với R Programming thì bạn đọc có thể viết một hàm có tên calculate_composite_index() để tính thống kê này với inputs đầu vào là các hệ số tải:

calculate_composite_index <- function(x) {
  
  ts <- sum(x)*sum(x)

  ms <- (ts + sum(1 - x^2))

  compositeReliability <- ts / ms
  
  return(compositeReliability)
  
}

Sử dụng hàm này chúng ta có thể tính toán Composite Reliability cho hai nhân tố DEL và WOM:

# For DEL: 
calculate_composite_index(all_loadings[6:8]) 
## [1] 0.8653525
# For WOM: 
calculate_composite_index(all_loadings[9:11]) 
## [1] 0.8637867

Đương nhiên vẫn có kết quả là 0.9486553 cho nhân tố ADA như cách tính tay:

calculate_composite_index(ada_loadings)
## [1] 0.9486553

Nếu không muốn viết hàm này thì có thể sử dụng bluegrafi để tính thống kê này.

Second Order CFA

Vẫn với bộ dữ liệu trên, xét mô hình Second Order CFA dưới đây:

Định nghĩa mô hình CFA này như sau:

my_SEM2 <- "EMP =~ empathy1 + empathy2 + empathy3 + empathy4 + empathy5
            SUR =~ surprise1 + surprise2
            WOM =~ WOM1 + WOM2 + WOM3
            EXP =~ Unique1 + Unique2 
            EXP ~ SUR 
            EXP ~ EMP
            SUR ~~ EMP
            WOM ~ EXP"

Chạy mô hình này đồng thời hiển thị các kết quả + các thống kê đánh giá độ phù hợp của mô hình:

# Run CFA model: 
SEM_model2 <- cfa(my_SEM2, data = data_f4.7)

# Show results: 
summary(SEM_model2, fit.measures = FALSE, standardized = FALSE)
## lavaan 0.6-9 ended normally after 39 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        28
##                                                       
##   Number of observations                           500
##                                                       
## Model Test User Model:
##                                                       
##   Test statistic                               146.058
##   Degrees of freedom                                50
##   P-value (Chi-square)                           0.000
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   EMP =~                                              
##     empathy1          1.000                           
##     empathy2          1.035    0.044   23.616    0.000
##     empathy3          0.949    0.041   23.240    0.000
##     empathy4          1.011    0.049   20.564    0.000
##     empathy5          0.917    0.047   19.439    0.000
##   SUR =~                                              
##     surprise1         1.000                           
##     surprise2         1.236    0.079   15.732    0.000
##   WOM =~                                              
##     WOM1              1.000                           
##     WOM2              1.075    0.042   25.722    0.000
##     WOM3              1.058    0.046   23.005    0.000
##   EXP =~                                              
##     Unique1           1.000                           
##     Unique2           0.977    0.070   13.925    0.000
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   EXP ~                                               
##     SUR               0.318    0.037    8.617    0.000
##     EMP               0.254    0.044    5.703    0.000
##   WOM ~                                               
##     EXP               0.866    0.095    9.084    0.000
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   EMP ~~                                              
##     SUR               0.301    0.041    7.419    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .empathy1          0.242    0.019   12.703    0.000
##    .empathy2          0.184    0.016   11.402    0.000
##    .empathy3          0.169    0.014   11.769    0.000
##    .empathy4          0.325    0.024   13.457    0.000
##    .empathy5          0.326    0.023   13.890    0.000
##    .surprise1         0.415    0.056    7.447    0.000
##    .surprise2         0.304    0.077    3.925    0.000
##    .WOM1              0.319    0.036    8.801    0.000
##    .WOM2              0.373    0.042    8.874    0.000
##    .WOM3              0.687    0.056   12.377    0.000
##    .Unique1           0.259    0.028    9.266    0.000
##    .Unique2           0.226    0.026    8.728    0.000
##     EMP               0.535    0.048   11.129    0.000
##     SUR               0.913    0.094    9.746    0.000
##    .WOM               0.938    0.083   11.344    0.000
##    .EXP               0.214    0.028    7.665    0.000
# Show some goodness of fit indices: 
fitMeasures(SEM_model2, 
            fit.measures = c("cfi", "tli", "rmsea", "gfi", "chisq", "df", "fmin", "ifi"))
##     cfi     tli   rmsea     gfi   chisq      df    fmin     ifi 
##   0.974   0.965   0.062   0.955 146.058  50.000   0.146   0.974

Model Diagram:

lavaanPlot(model = SEM_model2, coefs = TRUE, covs = TRUE, stars = c("regress"))

Final Notes

Post này hướng dẫn thực hiện CFA bậc 1 (First Order CFA) và bậc 2 (Second Order CFA) nhằm tại lập lại các kết quả của hai cases được trình bày trong textbook của J. Collier. Những gì được hướng dẫn trong post này là đã cover phần lớn những kết quả phổ biến thường được trình bày trong các peper sử dụng kĩ thuật phân tích CFA, SEM. Những kĩ thuật phân tích cũng như các vấn đề chưa đề cập khác nếu muốn bạn đọc có thể tham khảo từ cuốn sách trên của J. Collier.

Về công cụ sử dụng bạn đọc có thể dùng AMOS/Mplus (AMOS rất phổ biến ở VN) hoặc R (cụ thể là lavaan và các package hỗ trợ khác). Là người đã từng sử dụng những tools nêu trên, các nhân tôi thấy sử dụng R là thuận tiện hơn, và có thể thực hiện nhiều phân tích mà hiện giờ AMOS vẫn chưa có (dù sao thì các chức năng phân tích của AMOS hoặc các công cụ mất phí cùng loại thì các chức năng phân tích cũng là cố định và không thể tùy biến như R). Mặt khác cú pháp (syntax) rõ ràng và nhất quán của R/lavaan cũng làm cho việc thực hiện các phân tích CFA/SEM dễ thực hiện hơn.

References

  1. Structural Equation Modeling
  2. lavaan: An R Package for Structural Equation Modeling
  3. CONFIRMATORY FACTOR ANALYSIS (CFA) IN R WITH LAVAAN
  4. lavaan manual
  5. Applied Structural Equation Modeling using AMOS: Basic to Advanced Techniques
  6. Yang, Yanyun and Samuel B. Green. (2011), “Coefficient Alpha: A Reliability Coefficient for the 21st Century”, Journal of Psychoeducational Assessment, 29 (4), 377–392.
  7. Hair, Joseph F., William C. Black, Barry J. Babin, and Rolph E. Anderson. (2009), Multivariate Data Analysis (7th ed.). Upper Saddle River, NJ: Prentice Hall.