Phân tích dữ liệu bằng ngôn ngữ R - Viện Phương Nam (10-15/5/2025)
Thach Tran
2025-04-22
Ngày 6: Hồi qui logistic
Việc 1. Đọc dữ liệu ‘Bone data’
os = read.csv("C:\\Thach\\VN trips\\2025_2May\\Phuong Nam Institute\\Datasets\\Bone data.csv")
dim(os)## [1] 2162 9head(os)## id sex age weight height prior.fx fnbmd smoking fx
## 1 1 Male 73 98 175 0 1.08 1 0
## 2 2 Female 68 72 166 0 0.97 0 0
## 3 3 Male 68 87 184 0 1.01 0 0
## 4 4 Female 62 72 173 0 0.84 1 0
## 5 5 Male 61 72 173 0 0.81 1 0
## 6 6 Female 76 57 156 0 0.74 0 0Việc 2. Đánh giá mối liên quan giữa hút thuốc và nguy cơ gãy xương
m.1 = glm(fx ~ smoking, family = binomial, data = os)
summary(m.1)##
## Call:
## glm(formula = fx ~ smoking, family = binomial, data = os)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.03006 0.06442 -15.990 <2e-16 ***
## smoking -0.13796 0.10081 -1.368 0.171
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2441.4 on 2161 degrees of freedom
## Residual deviance: 2439.5 on 2160 degrees of freedom
## AIC: 2443.5
##
## Number of Fisher Scoring iterations: 4library(epiDisplay)## Warning: package 'epiDisplay' was built under R version 4.3.2## Loading required package: foreign## Loading required package: survival## Loading required package: MASS## Loading required package: nnetlogistic.display(m.1)##
## Logistic regression predicting fx
##
## OR(95%CI) P(Wald's test) P(LR-test)
## smoking: 1 vs 0 0.87 (0.71,1.06) 0.171 0.17
##
## Log-likelihood = -1219.7469
## No. of observations = 2162
## AIC value = 2443.4938[ChatGPT]
PROMPT: “Thực hiện mô hình hồi qui logistic đánh giá mối liên quan giữa hút thuốc và nguy cơ gãy xương.”
ChatGPT:
logistic_model <- glm(fx ~ smoking, data = df, family = binomial)
summary(logistic_model)
Việc 3. Đánh giá mối liên quan độc lập giữa hút thuốc và nguy cơ gãy xương
Giới tính là yếu tố gây nhiễu
library(compareGroups)
createTable(compareGroups(sex ~ smoking + fx, data = os))##
## --------Summary descriptives table by 'sex'---------
##
## _________________________________________
## Female Male p.overall
## N=1317 N=845
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## smoking 0.30 (0.46) 0.62 (0.49) <0.001
## fx 0.30 (0.46) 0.17 (0.38) <0.001
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Hiệu chỉnh nhiễu
os$gender = ifelse(os$sex == "Male", 1, 0)
m.2 = glm(fx ~ smoking + gender, family = binomial, data = os)
logistic.display(m.2)##
## Logistic regression predicting fx
##
## crude OR(95%CI) adj. OR(95%CI) P(Wald's test) P(LR-test)
## smoking: 1 vs 0 0.87 (0.71,1.06) 1.1 (0.89,1.36) 0.357 0.358
##
## gender: 1 vs 0 0.47 (0.38,0.58) 0.45 (0.36,0.57) < 0.001 < 0.001
##
## Log-likelihood = -1194.7997
## No. of observations = 2162
## AIC value = 2395.5994[ChatGPT]
PROMPT: “Y văn ghi nhận giới tính (sex) là yếu tố gây nhiễu. Bạn thực hiện mô hinh hồi qui logistic đa biến đánh giá mối liên quan độc lập giữa hút thuốc và nguy cơ gãy xương”
ChatGPT:
logistic_model_multivariate <- glm(fx ~ smoking + sex, data = df, family = binomial)
summary(logistic_model_multivariate)
Việc 4. Đánh giá tầm quan trọng của biến số
library(caret)## Loading required package: ggplot2##
## Attaching package: 'ggplot2'## The following object is masked from 'package:epiDisplay':
##
## alpha## Loading required package: lattice##
## Attaching package: 'lattice'## The following object is masked from 'package:epiDisplay':
##
## dotplot##
## Attaching package: 'caret'## The following object is masked from 'package:survival':
##
## clusterm.rel = glm(fx ~ gender + age + fnbmd + weight, family = binomial, data = os)
varImp(m.rel, scale = FALSE)## Overall
## gender 2.3022178
## age 2.5112611
## fnbmd 8.7193666
## weight 0.8239425[ChatGPT]
PROMPT 1: “Giả sử là mô hình dự báo nguy cơ gãy xương gồm tuổi, giới tính, mật độ xương ở cổ xương đùi và cân nặng. Hãy đánh giá tầm quan trọng tương đối (relative importance) của các biến số dự báo trên.”
ChatGPT: “Để đánh giá tầm quan trọng tương đối của các biến số dự báo trong mô hình hồi quy logistic, bạn có thể sử dụng phương pháp đánh giá tầm quan trọng của các biến số như phân tích tỷ lệ odds ratio, phương pháp bootstrap, hoặc mô hình cây quyết định (tree-based models) trong trường hợp bạn muốn tìm hiểu sự đóng góp của từng yếu tố.
Dưới đây là một cách để thực hiện việc này bằng cách sử dụng mô hình hồi quy logistic và phân tích tỷ lệ odds ratio (OR) của từng biến số. Odds ratio cho phép bạn đánh giá tầm quan trọng tương đối của mỗi yếu tố trong việc dự báo nguy cơ gãy xương.”
PROMPT 2: “Odds ratio không thể đánh giá mức độ quan trọng tương đối của biến số, Cụ thể với các biến định lượng (tuổi, mật độ xương và cân nặng) thì Odds ratio đánh giá thay đổi của odds gãy xương cho gia tăng mỗi đơn vị của biến định lượng; ngược lại với biến định tính như giới tính thì odds ratio đánh giá khác biệt về odds gãy xương giữa nữ so với nam. Bạn hãy đề nghị phương pháp thống kế khác tinh vi hơn”
ChatGPT: “Bạn hoàn toàn đúng! Odds Ratio không phải là phương pháp thích hợp nhất để đánh giá mức độ quan trọng tương đối của các biến số trong mô hình hồi quy logistic, đặc biệt đối với các biến định lượng như tuổi, mật độ xương và cân nặng. Để đánh giá tầm quan trọng của các biến này một cách chính xác hơn, có thể sử dụng các phương pháp khác như:
- **Phân tích tầm quan trọng tương đối bằng phương pháp randomForest hoặc caret Phương pháp này sử dụng mô hình cây quyết định ngẫu nhiên (Random Forest), giúp đánh giá tầm quan trọng của các biến dựa trên tác động của từng biến đến độ chính xác của mô hình khi biến đó bị loại bỏ.
-Random Forest tạo ra một bộ cây quyết định và sau đó đánh giá độ quan trọng của từng biến thông qua các kỹ thuật như mean decrease in accuracy (MDA) hoặc mean decrease in Gini (MDG).
Đây là phương pháp rất mạnh và linh hoạt, cho phép đánh giá tầm quan trọng của các biến số dựa trên dữ liệu thực nghiệm mà không yêu cầu giả định về phân phối của dữ liệu.
library(randomForest) library(caret)
Mô hình Random Forest để đánh giá tầm quan trọng các biến rf_model <- randomForest(fx ~ age + sex + fnbmd + weight, data = df)
Đánh giá tầm quan trọng của các biến importance(rf_model)
Vẽ biểu đồ tầm quan trọng của các biến varImpPlot(rf_model)
- Lasso (Least Absolute Shrinkage and Selection Operator) Phương pháp Lasso là một kỹ thuật hồi quy với phạt L1. Lasso giúp chọn lựa các biến quan trọng nhất trong mô hình hồi quy logistic bằng cách ép các hệ số của các biến không quan trọng về 0.
Với Lasso, bạn có thể chọn ra những biến có ảnh hưởng lớn nhất đến biến phụ thuộc (nguy cơ gãy xương trong trường hợp này).
Kỹ thuật này cực kỳ hữu ích khi có nhiều biến độc lập và bạn muốn xác định tầm quan trọng của từng biến trong việc dự báo nguy cơ.”
Việc 5. Xây dựng mô hình dự báo nguy cơ gãy xương
5.1 Tìm mô hình tối ưu bằng pp BMA
library(BMA)## Loading required package: leaps## Loading required package: robustbase##
## Attaching package: 'robustbase'## The following object is masked from 'package:survival':
##
## heart## Loading required package: inline## Loading required package: rrcov## Scalable Robust Estimators with High Breakdown Point (version 1.7-4)os.2 = na.omit(os)
xvars = os.2[, c("sex", "age", "weight", "height", "fnbmd", "smoking", "prior.fx")]
yvar = os.2[, c("fx")]
m.bma = bic.glm(xvars, yvar, strict = FALSE, OR = 20, glm.family = "binomial")
summary(m.bma)##
## Call:
## bic.glm.data.frame(x = xvars, y = yvar, glm.family = "binomial", strict = FALSE, OR = 20)
##
##
## 4 models were selected
## Best 4 models (cumulative posterior probability = 1 ):
##
## p!=0 EV SD model 1 model 2 model 3
## Intercept 100 2.249962 0.674258 2.509e+00 2.340e+00 1.140e+00
## sex 18.2 -0.048048 0.114551 . -2.578e-01 .
## age 16.2 0.002663 0.006838 . . 1.602e-02
## weight 0.0 0.000000 0.000000 . . .
## height 0.0 0.000000 0.000000 . . .
## fnbmd 100.0 -4.487636 0.436674 -4.595e+00 -4.279e+00 -4.302e+00
## smoking 0.0 0.000000 0.000000 . . .
## prior.fx 100.0 0.530729 0.134289 5.306e-01 5.445e-01 5.156e-01
##
## nVar 2 3 3
## BIC -1.402e+04 -1.402e+04 -1.402e+04
## post prob 0.696 0.142 0.122
## model 4
## Intercept 7.978e-01
## sex -2.851e-01
## age 1.780e-02
## weight .
## height .
## fnbmd -3.918e+00
## smoking .
## prior.fx 5.296e-01
##
## nVar 4
## BIC -1.402e+04
## post prob 0.040imageplot.bma(m.bma)
5.2 Mô hình tối ưu
os.2$fnbmd.sd = os.2$fnbmd/0.15
m.3 = glm(fx ~ fnbmd.sd, family = binomial, data = os.2)
logistic.display(m.3)##
## Logistic regression predicting fx
##
## OR(95%CI) P(Wald's test) P(LR-test)
## fnbmd.sd (cont. var.) 0.49 (0.43,0.54) < 0.001 < 0.001
##
## Log-likelihood = -1106.8757
## No. of observations = 2121
## AIC value = 2217.7514[ChatGPT]
PROMPT: “Xây dựng mô hình tối ưu dự báo gãy xương từ những biến số như tuổi (age), giới tính (sex), cân nặng (weight), chiều cao (height), mật độ xương (fnbmd), hút thuốc là (smoking) và tiền căn gãy xương (prior.fx). Dùng phương pháp Bayesian Model Averaging từ gói lệnh BMA”
ChatGPT:
Mô hình hồi quy logistic bma_model <- bic.glm(fx ~ age + sex + weight + height + fnbmd + smoking + prior.fx, data = df, glm.family = “binomial”)
In kết quả mô hình BMA summary(bma_model)