• Công thức kiểm định có dạng:

\[ HL = \sum_{i=1}^{k}\frac{n_i (p_i-\hat{p_i})^2}{\hat{p_i}(1-\hat{p_i})} \sim \chi^2(k-2) \] Kiểm định Hosmer-Lemeshow (1980) là kiểm định sự phù hợp của hàm hồi quy với biến độc lập là biến nhị phân. Kiểm định này được sử dụng rất phổ biến nhưng bên cạnh đó vẫn còn nhiều điểm chưa được phù hợp.

• Thông thường, để đánh giá mức độ giải thích (predict power) của biến độc lập cho biến phụ thuộc người ta dùng \(R^2\). Chỉ số này cho biết những biến độc lập giải thích được bao nhiêu % sự thay đổi của biến phụ thuộc.

• Ngược lại, goodness-of-fit (GOF) tests cho biết liệu mô hình đã chọn có chính xác hay không. Nếu p-value thấp (<0.05) loại bỏ mô hình, p-value cao thì chấp nhận mô hình.

• Mô hình với biến phụ thuộc là nhị phân, thông thường phải dùng link function (logit, probit, log-log hay hàm gì đó) tức là mô hình không trực tiếp dự báo xác suất. Vì vậy, kể cả trường hợp tuyến tính hay phi tuyến, vẫn có thể có trường hợp \(R^2\) cao, nhưng mô hình phân loại không chính xác và ngược lại \(R^2\) thấp nhưng mô hình lại phân loại chính xác

• Theo công thức HL mẫu được chia thành k nhóm bằng nhau, thông thường là 10 nhóm. Đối với mỗi nhóm, chúng ta tính tỷ lệ xảy ra sự kiện và kỳ vọng tỷ lệ xảy ra sự kiện, rồi lấy 2 kết quả trừ cho nhau.

• Ví dụ: sử dụng data “http://www.uam.es/personal_pdi/economicas/rsmanga/docs/mroz.dta”

mroz <- read_dta("D:/R/RMarkdown/hosmer_lemeshow/mroz_2.dta")

fit <- glm(inlf ~ kidslt6 + age + educ + huswage + exper, data = mroz, family = binomial(link = "logit") )

pred <- predict(fit, type = "response")

summary(fit)

## 
## Call:
## glm(formula = inlf ~ kidslt6 + age + educ + huswage + exper, 
##     family = binomial(link = "logit"), data = mroz)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5379  -0.9135   0.4501   0.8794   2.1043  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.25250    0.73787   1.697   0.0896 .  
## kidslt6     -1.45067    0.19887  -7.295 3.00e-13 ***
## age         -0.09754    0.01336  -7.301 2.86e-13 ***
## educ         0.21247    0.04226   5.027 4.98e-07 ***
## huswage     -0.04016    0.02113  -1.900   0.0574 .  
## exper        0.12120    0.01328   9.123  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1029.75  on 752  degrees of freedom
## Residual deviance:  815.22  on 747  degrees of freedom
## AIC: 827.22
## 
## Number of Fisher Scoring iterations: 4

f_hl <- function(gr){
  hoslem.test(mroz$inlf, fitted(fit), g = gr)$p.value
}


df_hl <- data.frame(gr = 5:25, hl = sapply(5:25, f_hl)) %>% 
  mutate(reject = case_when(hl < 0.05 ~ "Lower", TRUE ~ "Upper") %>% as.factor())

df_hl %>% 
ggplot(aes(x = gr, y = hl, color = reject))+
  geom_point(size = 4)+
  labs(title = "Correlation between group vs p-value of HL test",
       subtitle = "S. Lemeshow and D.W. Hosmer (1982)",
       y = "p-value",
       x = "No of group")+
  geom_hline(yintercept = 0.05, color = "purple")+
  scale_y_continuous(breaks = c(0, 0.05, 0.2, 0.4, 0.6, 0.8))+
  theme(legend.position = "none")

Đường gạch ngang là mức xác suất p-value = 0.05, nếu dựa vào mức xác suất này để kết luận thì sẽ không chính xác. Vì khi thay đổi số nhóm thì p-value cũng thay đổi như trên hình vẽ.

• Năm 1997, S.le Cessie cũng có đề xuất 1 chút thay đổi trong công thức kiểm định. Công thức này dựa trên ý tưởng làm trơn phần dư (smoothed residual based).

  • Ưu điểm của công thức HL (1982) là ước lượng p-value dựa trên việc phân chia thành các nhóm cố định. Vì vậy, nó rất dễ hiểu. Nhưng nhược điểm là p-value phụ thuộc vào việc lựa chọn số lượng nhóm và điểm cutpoints.

  • le Cessie and van Houwelingen đề xuất việc phân nhóm dựa trên smoothed residuals và sử dụng weight function đối với từng nhóm

Kết quả của kiểm định này được thể hiện trong đồ thị dưới

# Hàm với 
library(MKmisc)
f_hl_2 <- function(gr){
  HLgof.test(fit = fitted(fit), obs = mroz$inlf, ngr = gr)$H$p.value 
}

df_hl <- data.frame(gr = 5:25, hl = sapply(5:25, f_hl_2)) %>% 
  mutate(reject = case_when(hl < 0.05 ~ "Lower", TRUE ~ "Upper") %>% as.factor())

df_hl %>% 
ggplot(aes(x = gr, y = hl, color = reject))+
  geom_point(size = 4)+
  labs(title = "Correlation between group vs p-value of HL test",
       subtitle = "D.W. Hosmer, T. Hosmer, S. le Cessie, S. Lemeshow (1997)",
       y = "p-value",
       x = "No of group")+
  geom_hline(yintercept = 0.05, color = "purple")+
  scale_y_continuous(breaks = c(0, 0.05, 0.2, 0.4, 0.6, 0.8))+
  scale_color_brewer(palette = "Set1")+
  theme(legend.position = "none")

Kết quả dường như cũng không tốt hơn nhiều so với hàm HL (1980)