Giới thiệu
Bài này chúng ta sẽ tái tạo lại hầu hết các kết quả cho chương 5 cuốn Marketing Research with SPSS. Những tính toán thủ công này giúp chúng ta hiểu rõ hơn cách thức tính toán - phân tích dữ liệu. Ngoài ra, chúng ta cũng thấy rằng SPSS có sai sót.
Điều này không có gì lạ: SPSS không phải là duy nhất có tình trạng này. Ví dụ, Stata cho đến phiên bản 14 vẫn báo cáo sai R2 khi phân tích dữ liệu mảng - điều mà giới thống kê và phân tích dữ liệu đã chỉ ra nhưng công ti sản xuất ra Stata vẫn chưa chỉnh sửa dù lỗi này được phát hiện từ lâu.
# Đọc dữ liệu
path <- dir("F:/marketing", full.names = TRUE)
path## [1] "F:/marketing/ch01" "F:/marketing/ch02" "F:/marketing/ch03"
## [4] "F:/marketing/ch04" "F:/marketing/ch05" "F:/marketing/ch06"
## [7] "F:/marketing/ch07" "F:/marketing/ch08" "F:/marketing/ch09"
## [10] "F:/marketing/ch10" "F:/marketing/ch11"du_lieu <- dir("F:/marketing/ch05", full.names = TRUE)
du_lieu## [1] "F:/marketing/ch05/pizza.sav"
## [2] "F:/marketing/ch05/pizza_incl_gender.sav"
## [3] "F:/marketing/ch05/pizza_stepwise_method.sav"library(foreign)
org_df <- read.spss(du_lieu[1], to.data.frame = TRUE)
# Xem qua dữ liệu:
library(tidyverse)
org_df %>% head()## resp_nr satisfaction reception service
## 1 16 very dissatisfied dissatisfied very satisfied
## 2 14 dissatisfied dissatisfied neutral
## 3 28 dissatisfied dissatisfied rather dissatisfied
## 4 58 dissatisfied dissatisfied rather dissatisfied
## 5 7 rather dissatisfied neutral rather satisfied
## 6 17 rather dissatisfied neutral rather dissatisfied
## waitingtime quality price
## 1 very dissatisfied dissatisfied dissatisfied
## 2 very dissatisfied very dissatisfied very dissatisfied
## 3 rather dissatisfied dissatisfied very dissatisfied
## 4 rather dissatisfied dissatisfied very dissatisfied
## 5 satisfied rather dissatisfied rather satisfied
## 6 rather dissatisfied dissatisfied very dissatisfiedorg_df %>% dim()## [1] 107 7org_df %>% str()## 'data.frame': 107 obs. of 7 variables:
## $ resp_nr : num 16 14 28 58 7 17 29 35 36 50 ...
## $ satisfaction: Factor w/ 7 levels "very dissatisfied",..: 1 2 2 2 3 3 3 3 3 3 ...
## $ reception : Factor w/ 7 levels "very dissatisfied",..: 2 2 2 2 4 4 3 4 3 3 ...
## $ service : Factor w/ 7 levels "very dissatisfied",..: 7 4 3 3 5 3 5 5 4 6 ...
## $ waitingtime : Factor w/ 7 levels "very dissatisfied",..: 1 1 3 3 6 3 5 4 4 7 ...
## $ quality : Factor w/ 7 levels "very dissatisfied",..: 2 1 2 2 3 2 4 3 3 3 ...
## $ price : Factor w/ 7 levels "very dissatisfied",..: 2 1 1 1 5 1 2 1 2 1 ...
## - attr(*, "variable.labels")= Named chr "Respondent number" "Overall satisfaction" "Reception" "Service" ...
## ..- attr(*, "names")= chr "resp_nr" "satisfaction" "reception" "service" ...# Chuyển hóa các biến factor về numeric:
df1 <- org_df %>% mutate_if(is.factor, as.numeric)
# loại ra biến không cần thiết:
df1 <- df1 %>% select(-1)
#----------------------------------------
# Thực hiện hồi quy OLS không chuẩn hóa:
#----------------------------------------
ols1 <- lm(satisfaction ~ ., data = df1)
# Tái lập một phần kết quả ở Figure 5.32:
ols1 %>% summary()##
## Call:
## lm(formula = satisfaction ~ ., data = df1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4639 -0.4693 -0.0590 0.6420 1.6930
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.17164 0.54752 -0.313 0.75456
## reception 0.23080 0.08979 2.571 0.01161 *
## service 0.26477 0.09878 2.680 0.00859 **
## waitingtime 0.17947 0.08520 2.106 0.03765 *
## quality 0.44151 0.10561 4.181 6.2e-05 ***
## price 0.12903 0.05930 2.176 0.03190 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8589 on 101 degrees of freedom
## Multiple R-squared: 0.5753, Adjusted R-squared: 0.5543
## F-statistic: 27.36 on 5 and 101 DF, p-value: < 2.2e-16# Tính các VIF như Figure 3.32:
library(car)
vif(ols1)## reception service waitingtime quality price
## 1.365974 1.391428 1.649332 1.333917 1.552188#-----------------------------------------
# Viết hàm chuẩn hóa
#-----------------------------------------
# Do chúng ta có nhiều phân tích dính đến chữ
# "chuẩn hóa" nên chúng ta viết một hàm chuẩn hóa
# để thực hiện hồi quy chuẩn hóa nếu cần:
chuan_hoa <- function(x) {
ch <- (x - mean(x)) / sd(x)
return(ch)
}
# Ngoài ra còn viết hàm chuẩn hóa phần dư theo
# công thức http://www.r-tutor.com/elementary-statistics/simple-linear-regression/standardized-residual
# như sau:
chuan_hoa_phan_du <- function(y) {
ch_phan_du <- y / sd(y)
return(ch_phan_du)
}
# Tái tạo lại hình ảnh như Figure 5.21. Trước
# hết lấy ra phần dư:
phan_du <- ols1$residuals
# Rồi chuẩn hóa phần dư:
phan_du_ch <- chuan_hoa_phan_du(phan_du)
# Chú ý hai đặc điểm sau của các biến được chuẩn hóa:
# (1) mean luôn là 0 , và (2) độ lệch chuẩn luôn là 1.
# Có thể kiếm tra lại đặc điểm này:
mean(phan_du_ch)## [1] -1.950956e-17sd(phan_du_ch)## [1] 1# Có thể tạo ra Histogram như Figure 5.21. Trước hết
# biến phần dư được chuẩn hóa về data frame:
df_phan_du <- as.data.frame(phan_du_ch)
# Rồi vẽ. Chú ý rằng do số bin là khác nên hình vẽ
# cho Histogram không nhất thiết phải giống. Cái mà
# chúng ta cần là từ hình ảnh của chúng để chẩn
# đoán phần dư chuẩn hóa (và cả phần dư không chuẩn hóa)
# có phân phối chuẩn hay không:
theme_set(theme_minimal())
df_phan_du %>%
ggplot(aes(phan_du_ch)) +
geom_density(color = "red", fill = "red", alpha = 0.3) +
geom_histogram(aes(y = ..density..), color = "blue", fill = "blue", alpha = 0.3)
# Vầ QQ Plot như Figure 5.21. Hình này thì không có gì khác biệt:
df_phan_du %>%
ggplot(aes(sample = phan_du_ch)) +
geom_abline(intercept = 0, slope = 1, color = "red") +
stat_qq()
# Ý nghĩa của hai hình này là: chúng đều chỉ ra bằng
# chứng mạnh mẽ rằng phần dư chuẩn hóa (và cả không chuẩn)
# là phân phối chuẩn. Đây là một giả thuyết quan trọng
# của mô hình hồi quy OLS cổ điển. Chúng ta có thể đưa
# ra bằng chứng thống kê rằng phần dư là phân phối chuẩn
# bằng Shapiro - Wilk để tái tạo lại như ở Figure 2.25:
shapiro.test(phan_du)##
## Shapiro-Wilk normality test
##
## data: phan_du
## W = 0.98556, p-value = 0.3024#------------------------------------
# Thực hiện hồi quy chuẩn hóa
#------------------------------------
# Trước hết chuẩn hóa dữ liệu:
df1_ch <- df1 %>% mutate_all(chuan_hoa)
# Rồi thực hiện hồi quy chuẩn hóa:
ols_ch <- lm(satisfaction ~ ., data = df1_ch)
# Tái tạo lại kết quả ở Figure 5.32:
ols_ch %>% summary()##
## Call:
## lm(formula = satisfaction ~ ., data = df1_ch)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.91529 -0.36478 -0.04586 0.49905 1.31601
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.292e-16 6.454e-02 0.000 1.00000
## reception 1.948e-01 7.579e-02 2.571 0.01161 *
## service 2.050e-01 7.649e-02 2.680 0.00859 **
## waitingtime 1.754e-01 8.328e-02 2.106 0.03765 *
## quality 3.131e-01 7.490e-02 4.181 6.2e-05 ***
## price 1.758e-01 8.079e-02 2.176 0.03190 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6676 on 101 degrees of freedom
## Multiple R-squared: 0.5753, Adjusted R-squared: 0.5543
## F-statistic: 27.36 on 5 and 101 DF, p-value: < 2.2e-16# Thực tế là chúng ta đang thực hiện hồi quy không
# có hệ số chặn nên có thể làm như sau:
ols_ch_cach_khac <- lm(satisfaction ~ 0 + ., data = df1_ch)
# Kết quả là không có gì khác biệt:
ols_ch_cach_khac %>% summary()##
## Call:
## lm(formula = satisfaction ~ 0 + ., data = df1_ch)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.91529 -0.36478 -0.04586 0.49905 1.31601
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## reception 0.19482 0.07542 2.583 0.01121 *
## service 0.20503 0.07612 2.694 0.00826 **
## waitingtime 0.17542 0.08287 2.117 0.03671 *
## quality 0.31311 0.07453 4.201 5.7e-05 ***
## price 0.17579 0.08039 2.187 0.03106 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6644 on 102 degrees of freedom
## Multiple R-squared: 0.5753, Adjusted R-squared: 0.5545
## F-statistic: 27.63 on 5 and 102 DF, p-value: < 2.2e-16# Tạo ra hình như Figure 5.26. Trước hết lấy ra
# các giá trị dự báo:
du_bao <- ols1$fitted.values
# Rồi chuẩn hóa nó:
du_bao_ch <- du_bao %>% chuan_hoa()
df2 <- data.frame(phan_du_ch, du_bao_ch)
df2 %>% ggplot(aes(du_bao_ch, phan_du_ch)) + geom_point()
# Các hình từ Figure 5.26 đến 5.27 được tạo ra với mục đích
# duy nhất là kiểm tra về mặt hình ảnh của 1 trong các giả
# thuyết của mô hình OLS tuyến tính cổ điển là điều kiện
# phương sai không thay đổi (homoscedaticity). Các hình này
# là nhiều một cách không cần thiết. Để kiểm tra điều này
# chỉ cần 1 plot mà thôi (xem 9.2.1.2 tài liệu KTL ứng dụng với R).
# Chẳng hạn chúng ta có thể plot phần dư với, ví dụ, biến service:
df3 <- df1 %>% mutate(Resid = phan_du^2)
df3 %>% ggplot(aes(service, Resid)) + geom_point()
# Hoặc với giá trị dự báo:
df3 <- df3 %>% mutate(Pred2 = ols1$fitted.values)
df3 %>% ggplot(aes(Pred2, Resid)) + geom_point()
# Nhận xét: Các hình ảnh này cho thấy có thể giả thuyết
# phương sai đồng nhất là bị vi phạm. Kiểm định chính thức
# các bạn có thể xem ở chương 9 tài liệu KTL ứng dụng với R.
# Sử dụng hình ảnh để chẩn đoán các giả thiết bị vi phạm
# chúng ta thực tế chỉ cần một lệnh ngắn gọn mà thôi:
library("ggfortify")
autoplot(ols1, which = 1:6, ncol = 3, label.size = 2)
# Hình ảnh này cho thấy, ví dụ, nếu căn cứ tiêu chí
# Cook's Distance thì quan sát thứ 37 và 62 là các outliers.
# Tuy nhiên trong sách này thì người ta lại sử dụng một
# cách tiếp cận khác: bất cứ quan sát nào có phần dư tương
# ứng nằm ngoài 2 lần độ lệch chuẩn từ mean của phần dư
# thì là outliers. Tính toán của SPSS chỉ ra đó là các quan
# sát thứ 1, 5, 7, 10, và 16. Chúng ta có thể tái lập lại
# những kết quả này như sau:
ct <- mean(phan_du) + 2*sd(phan_du) # Cận trên của 2 lần lệch chuẩn từ mean.
cd <- mean(phan_du) - 2*sd(phan_du) # Cận dưới.
# Ghép thêm phần dư vào:
df1 <- df1 %>% mutate(PhanDu = phan_du)
# Chỉ ra những quan sát nào là outliers theo tiêu chuẩn trên:
df1 <- df1 %>% mutate(Outlies = case_when(PhanDu > ct | PhanDu < cd ~ "Yes",
PhanDu <= ct | PhanDu >= cd ~ "No"))
# Tạo thêm ID để thấy rằng thứ mà chúng ta hiểu và tính toán
# để so lại với SPSS:
df1 <- df1 %>% mutate(ID = 1:nrow(df1))
# Các phần tử bất thường này là:
df1 %>% filter(Outlies == "Yes")## satisfaction reception service waitingtime quality price PhanDu
## 1 1 2 7 1 2 2 -2.463896
## 2 3 4 5 6 3 5 -2.121921
## 3 3 3 5 5 4 2 -1.766058
## 4 3 3 6 7 3 1 -1.819229
## 5 4 5 6 7 3 5 -1.796961
## 6 5 3 3 3 3 1 1.692962
## Outlies ID
## 1 Yes 1
## 2 Yes 5
## 3 Yes 7
## 4 Yes 10
## 5 Yes 16
## 6 Yes 37# Kết quả tính toán của chúng ta là gần như trùng khớp
# với cách chơi của SPSS. Câu hỏi là ai đúng? Câu trả lời
# là R đúng chứ không phải SPSS. Có hai bằng chứng để
# tự tin vào điều này. Thứ nhất, quan sát thứ 37 có trong
# danh sách đen từ plot ơ trên. Thứ hai (và cũng là chắc)
# chắn nhất: phần dư của phần tử thứ 37 là 1.692962 lớn hơn
# cận trên 1.676749 như chúng ta có thể thấy. Điều này có
# nghĩa là SPSS bỏ sót chính một phần tử theo tiêu chuẩn
# nằm ngoài 2 lần độ lệch chuẩn từ mean phần dư mà nó định nghĩa:
ct## [1] 1.676749(Còn nữa… )