#Ngày 6: Hồi qui logistic
#Việc 1. Dữ liệu "Bone data.csv" và gọi dữ liệu này là “os”
dbone <- read.csv("E:\\HOPT\\PTDLR\\DATA\\Bone data.csv")
head(dbone)
## 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 0
dim(dbone)
## [1] 2162 9
#Việc 2. Đánh giá mối liên quan giữa hút thuốc (smoking) và nguy cơ gãy xương (fx)
#2.1 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. Viết 1 câu văn giải thích kết quả
m.1 = glm(fx ~ smoking, family = binomial, data = dbone)
summary(m.1)
##
## Call:
## glm(formula = fx ~ smoking, family = binomial, data = dbone)
##
## 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: 4
#thư viện epiDisplay dành cho phân tích dữ liệu y tế và dịch tễ học (epidemiology). Nó cung cấp các hàm thân thiện,
#Mục đích Ví dụ hàm
#1. Tóm tắt biến định tính (categorical) tab1()
#2.Tóm tắt biến định lượng (numeric) tabstat()
#3. Thống kê mô tả toàn bộ dữ liệu des()
#4.Phân tích hồi quy logistic logistic.display()
#5.Phân tích hồi quy tuyến tính regress.display()
#6Kiểm định thống kê đơn giản ci.mean() – tính CI, ci.prop()
library(epiDisplay)
## Loading required package: foreign
## Loading required package: survival
## Loading required package: MASS
## Loading required package: nnet
logistic.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
#Kiểm tra biến fx (gãy xương) và smoke (hút thuốc). Nếu cần, chuyển về factor:
#dbone$fx <- factor(dbone$fx, levels = c(0, 1), labels = c("Không gãy", "Gãy"))
#dbone$smoke <- factor(dbone$smoking, levels = c(0, 1), labels = c("Không hút", "Hút thuốc"))
# 1. Đọc dữ liệu
dbone <- read.csv("E:/HOPT/PTDLR/DATA/Bone data.csv") # chỉnh đường dẫn nếu cần
# 2. Kiểm tra tên và định dạng biến
str(dbone)
## 'data.frame': 2162 obs. of 9 variables:
## $ id : int 1 2 3 4 5 6 7 8 9 10 ...
## $ sex : chr "Male" "Female" "Male" "Female" ...
## $ age : int 73 68 68 62 61 76 63 64 76 64 ...
## $ weight : int 98 72 87 72 72 57 97 85 48 89 ...
## $ height : int 175 166 184 173 173 156 173 167 153 166 ...
## $ prior.fx: int 0 0 0 0 0 0 0 0 0 0 ...
## $ fnbmd : num 1.08 0.97 1.01 0.84 0.81 0.74 1.01 0.86 0.65 0.98 ...
## $ smoking : int 1 0 0 1 1 0 1 0 0 0 ...
## $ fx : int 0 0 0 0 0 0 0 0 0 0 ...
names(dbone)
## [1] "id" "sex" "age" "weight" "height" "prior.fx" "fnbmd"
## [8] "smoking" "fx"
# 3. Kiểm tra phân bố biến smoking
table(dbone$smoking, useNA = "always")
##
## 0 1 <NA>
## 1243 919 0
# 4. Tạo biến phân loại từ biến smoking
# Nếu smoking là 0/1:
dbone$smoking <- factor(dbone$smoking, levels = c(0, 1),labels = c("Không hút", "Hút thuốc"))
# 5. Kiểm tra lại sau khi chuyển đổi
table(dbone$smoking)
##
## Không hút Hút thuốc
## 1243 919
# 6. Chuyển fx thành biến nhị phân dạng factor nếu chưa
dbone$fx <- factor(dbone$fx, levels = c(0, 1),labels = c("Không gãy", "Gãy"))
# 7. Xây dựng mô hình hồi quy logistic
model_logit <- glm(fx ~ smoking, data = dbone, family = binomial)
# 8. Xem kết quả
summary(model_logit)
##
## Call:
## glm(formula = fx ~ smoking, family = binomial, data = dbone)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.03006 0.06442 -15.990 <2e-16 ***
## smokingHút thuốc -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: 4
# 9. (Tùy chọn) Hiển thị mô hình đẹp hơn bằng gói epiDisplay
# install.packages("epiDisplay") # nếu chưa có
library(epiDisplay)
logistic.display(model_logit)
##
## Logistic regression predicting fx : Gãy vs Không gãy
##
## OR(95%CI) P(Wald's test) P(LR-test)
## smoking: Hút thuốc vs Không hút 0.87 (0.71,1.06) 0.171 0.17
##
## Log-likelihood = -1219.7469
## No. of observations = 2162
## AIC value = 2443.4938
#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
logistic_model <- glm(fx ~ smoking, data = dbone, family = binomial)
summary(logistic_model)
##
## Call:
## glm(formula = fx ~ smoking, family = binomial, data = dbone)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.03006 0.06442 -15.990 <2e-16 ***
## smokingHút thuốc -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: 4
#Việc 3. Đánh giá mối liên quan độc lập giữa hút thuốc (smoking) và nguy cơ gãy xương (fx)sau khi hiệu chỉnh cho ảnh hưởng nhiễu của giới tính (sex)
#3.1 Thực hiện mô hình hối qui logistic đá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 = dbone))
##
## --------Summary descriptives table by 'sex'---------
##
## _______________________________________________
## Female Male p.overall
## N=1317 N=845
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## smoking: <0.001
## Không hút 923 (70.1%) 320 (37.9%)
## Hút thuốc 394 (29.9%) 525 (62.1%)
## fx: <0.001
## Không gãy 916 (69.6%) 701 (83.0%)
## Gãy 401 (30.4%) 144 (17.0%)
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
#Hiệu chỉnh nhiễu
dbone$gender = ifelse(dbone$sex == "Male", 1, 0)
m.2 = glm(fx ~ smoking + gender, family = binomial, data = dbone)
logistic.display(m.2)
##
## Logistic regression predicting fx : Gãy vs Không gãy
##
## crude OR(95%CI) adj. OR(95%CI)
## smoking: Hút thuốc vs Không hút 0.87 (0.71,1.06) 1.1 (0.89,1.36)
##
## gender: 1 vs 0 0.47 (0.38,0.58) 0.45 (0.36,0.57)
##
## P(Wald's test) P(LR-test)
## smoking: Hút thuốc vs Không hút 0.357 0.358
##
## gender: 1 vs 0 < 0.001 < 0.001
##
## Log-likelihood = -1194.7997
## No. of observations = 2162
## AIC value = 2395.5994
#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
logistic_model_multivariate <- glm(fx ~ smoking + sex, data = dbone, family = binomial)
summary(logistic_model_multivariate)
##
## Call:
## glm(formula = fx ~ smoking + sex, family = binomial, data = dbone)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.85598 0.06835 -12.523 < 2e-16 ***
## smokingHút thuốc 0.09872 0.10724 0.921 0.357
## sexMale -0.78880 0.11493 -6.863 6.73e-12 ***
## ---
## 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: 2389.6 on 2159 degrees of freedom
## AIC: 2395.6
##
## Number of Fisher Scoring iterations: 4
#Việc 4. 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.
#4.1 Hãy đánh giá tầm quan trọng của các biến số dự báo trên.
#cài caret mục đích
# Huấn luyện mô hình dự báo Classification (phân loại) và Regression (hồi quy)
#Hỗ trợ hơn 200 mô hình Hỗ trợ từ glm, rpart, svm, rf, xgb, v.v.
# Tiền xử lý dữ liệu Chuẩn hóa, xử lý NA, chuyển đổi biến, PCA
# Resampling tự động Cross-validation, bootstrap, LOOCV
# So sánh mô hình Tự động đánh giá độ chính xác các thuật toán
#install.packages("caret")
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':
##
## cluster
#xử lý nhiểu data
str(dbone$fx)
## Factor w/ 2 levels "Không gãy","Gãy": 1 1 1 1 1 1 1 1 1 1 ...
table(dbone$fx)
##
## Không gãy Gãy
## 1617 545
# Chuyển fx về kiểu nhị phân numeric
# Tạo lại biến nhị phân
dbone$fx_bin <- ifelse(dbone$fx == "Gãy", 1,
ifelse(dbone$fx == "Không gãy", 0, NA))
dbone$fx_bin <- as.numeric(as.character(dbone$fx))
## Warning: NAs introduced by coercion
dbone_clean <- na.omit(dbone[, c("fx_bin", "gender", "age", "fnbmd", "weight")])
table(dbone$fx_bin, useNA = "always")
##
## <NA>
## 2162
nrow(dbone_clean) # Phải > 0
## [1] 0
colSums(is.na(dbone[, c("fx_bin", "gender", "age", "fnbmd", "weight")]))
## fx_bin gender age fnbmd weight
## 2162 0 0 40 1
table(dbone$fx, useNA = "always")
##
## Không gãy Gãy <NA>
## 1617 545 0
str(dbone$fx)
## Factor w/ 2 levels "Không gãy","Gãy": 1 1 1 1 1 1 1 1 1 1 ...
#fix all
# Đọc dữ liệu
# Với Windows:
dbone <- read.csv("E:\\HOPT\\PTDLR\\DATA\\Bone data.csv", fileEncoding = "UTF-8-BOM")
# Kiểm tra biến fx
table(dbone$fx, useNA = "always")
##
## 0 1 <NA>
## 1617 545 0
str(dbone$fx)
## int [1:2162] 0 0 0 0 0 0 0 0 0 0 ...
dbone$gender <- ifelse(dbone$sex == "Male", 1, 0)
dbone_clean <- na.omit(dbone[, c("fx", "gender", "age", "fnbmd", "weight")])
m.rel <- glm(fx ~ gender + age + fnbmd + weight, family = binomial, data = dbone_clean)
summary(m.rel)
##
## Call:
## glm(formula = fx ~ gender + age + fnbmd + weight, family = binomial,
## data = dbone_clean)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.806796 0.770234 1.047 0.2949
## gender -0.294860 0.128077 -2.302 0.0213 *
## age 0.019647 0.007824 2.511 0.0120 *
## fnbmd -4.322857 0.495777 -8.719 <2e-16 ***
## weight 0.004078 0.004949 0.824 0.4100
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2395.8 on 2120 degrees of freedom
## Residual deviance: 2203.3 on 2116 degrees of freedom
## AIC: 2213.3
##
## Number of Fisher Scoring iterations: 4
library(caret)
m.rel = glm(fx ~ gender + age + fnbmd + weight, family = binomial, data = dbone)
varImp(m.rel, scale = FALSE)
## Overall
## gender 2.3022178
## age 2.5112611
## fnbmd 8.7193666
## weight 0.8239425
#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 dự báo nguy cơ gãy xương bằng phương pháp 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-7)
dbone2 = na.omit(dbone)
xvars = dbone2[, c("sex", "age", "weight", "height", "fnbmd", "smoking", "prior.fx")]
yvar = dbone2[, 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.040
imageplot.bma(m.bma)
#5.2 Viết phương trình của mô hình dự báo
install.packages("epiDisplay") # chỉ cần 1 lần
## Warning: package 'epiDisplay' is in use and will not be installed
library(epiDisplay)
dbone2$fnbmd.sd = dbone2$fnbmd/0.15
m.3 = glm(fx ~ fnbmd.sd, family = binomial, data = dbone2)
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
#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”
# Gọi thư viện nếu chưa
library(BMA)
# Mô hình BMA
bma_model <- bic.glm(fx ~ age + sex + weight + height + fnbmd + smoking + prior.fx,
data = dbone,
glm.family = "binomial")
# Hiển thị kết quả
summary(bma_model)
##
## Call:
## bic.glm.formula(f = fx ~ age + sex + weight + height + fnbmd + smoking + prior.fx, data = dbone, glm.family = "binomial")
##
##
## 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
## age 16.2 0.002663 0.006838 . . 1.602e-02
## sexMale 18.2 -0.048048 0.114551 . -2.578e-01 .
## 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
## age 1.780e-02
## sexMale -2.851e-01
## weight .
## height .
## fnbmd -3.918e+00
## smoking .
## prior.fx 5.296e-01
##
## nVar 4
## BIC -1.402e+04
## post prob 0.040
##
## 41 observations deleted due to missingness.