1. Hồi qui logistic – hiệu chỉnh yếu tố nhiễu

1.1 Đọc dữ liệu

df <- read.csv("D:/RData/Stroke Data.csv")

dim(df)

## [1] 5110   12

1.2 Mô tả dữ liệu theo tình trạng đột quỵ

library(table1)

## 
## Attaching package: 'table1'

## The following objects are masked from 'package:base':
## 
##     units, units<-

table1(~ . | stroke, data = df)

## Warning in table1.formula(~. | stroke, data = df): Terms to the right of '|' in
## formula 'x' define table columns and are expected to be factors with meaningful
## labels.

	0 (N=4861)	1 (N=249)	Overall (N=5110)
id
Mean (SD)	36500 (21100)	37100 (22000)	36500 (21200)
Median [Min, Max]	37000 [67.0, 72900]	36700 [210, 72900]	36900 [67.0, 72900]
gender
Female	2853 (58.7%)	141 (56.6%)	2994 (58.6%)
Male	2008 (41.3%)	108 (43.4%)	2116 (41.4%)
age
Mean (SD)	42.0 (22.3)	67.7 (12.7)	43.2 (22.6)
Median [Min, Max]	43.0 [0.0800, 82.0]	71.0 [1.32, 82.0]	45.0 [0.0800, 82.0]
hypertension
Mean (SD)	0.0889 (0.285)	0.265 (0.442)	0.0975 (0.297)
Median [Min, Max]	0 [0, 1.00]	0 [0, 1.00]	0 [0, 1.00]
heart.disease
Mean (SD)	0.0471 (0.212)	0.189 (0.392)	0.0540 (0.226)
Median [Min, Max]	0 [0, 1.00]	0 [0, 1.00]	0 [0, 1.00]
ever.married
No	1728 (35.5%)	29 (11.6%)	1757 (34.4%)
Yes	3133 (64.5%)	220 (88.4%)	3353 (65.6%)
work.type
children	685 (14.1%)	2 (0.8%)	687 (13.4%)
Govt_job	624 (12.8%)	33 (13.3%)	657 (12.9%)
Never_worked	22 (0.5%)	0 (0%)	22 (0.4%)
Private	2776 (57.1%)	149 (59.8%)	2925 (57.2%)
Self-employed	754 (15.5%)	65 (26.1%)	819 (16.0%)
Residence.type
Rural	2400 (49.4%)	114 (45.8%)	2514 (49.2%)
Urban	2461 (50.6%)	135 (54.2%)	2596 (50.8%)
glucose.level
Mean (SD)	105 (43.8)	133 (61.9)	106 (45.3)
Median [Min, Max]	91.5 [55.1, 268]	105 [56.1, 272]	91.9 [55.1, 272]
bmi
Mean (SD)	28.8 (7.91)	30.5 (6.33)	28.9 (7.85)
Median [Min, Max]	28.0 [10.3, 97.6]	29.7 [16.9, 56.6]	28.1 [10.3, 97.6]
Missing	161 (3.3%)	40 (16.1%)	201 (3.9%)
smoking
formerly smoked	815 (16.8%)	70 (28.1%)	885 (17.3%)
never smoked	1802 (37.1%)	90 (36.1%)	1892 (37.0%)
smokes	747 (15.4%)	42 (16.9%)	789 (15.4%)
Unknown	1497 (30.8%)	47 (18.9%)	1544 (30.2%)

1.3 Mối liên quan giữa cao huyết áp và đột quỵ

str(df)

## 'data.frame':    5110 obs. of  12 variables:
##  $ id            : int  67 77 84 91 99 121 129 132 156 163 ...
##  $ gender        : chr  "Female" "Female" "Male" "Female" ...
##  $ age           : num  17 13 55 42 31 38 24 80 33 20 ...
##  $ hypertension  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ heart.disease : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ ever.married  : chr  "No" "No" "Yes" "No" ...
##  $ work.type     : chr  "Private" "children" "Private" "Private" ...
##  $ Residence.type: chr  "Urban" "Rural" "Urban" "Urban" ...
##  $ glucose.level : num  93 85.8 89.2 98.5 108.9 ...
##  $ bmi           : num  NA 18.6 31.5 18.5 52.3 NA 26.2 NA 42.2 28.8 ...
##  $ smoking       : chr  "formerly smoked" "Unknown" "never smoked" "never smoked" ...
##  $ stroke        : int  0 0 0 0 0 0 0 0 0 0 ...

modelLog_Hypertension <- glm(
  stroke ~ hypertension,
  family = binomial,
  data = df
)

summary(modelLog_Hypertension)

## 
## Call:
## glm(formula = stroke ~ hypertension, family = binomial, data = df)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.18644    0.07543 -42.242   <2e-16 ***
## hypertension  1.30767    0.15217   8.593   <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: 1990.4  on 5109  degrees of freedom
## Residual deviance: 1929.3  on 5108  degrees of freedom
## AIC: 1933.3
## 
## Number of Fisher Scoring iterations: 6

coef(modelLog_Hypertension)

##  (Intercept) hypertension 
##    -3.186443     1.307672

exp(coef(modelLog_Hypertension))

##  (Intercept) hypertension 
##   0.04131858   3.69755616

exp(confint(modelLog_Hypertension))

## Waiting for profiling to be done...

##                   2.5 %     97.5 %
## (Intercept)  0.03551723 0.04774556
## hypertension 2.72798259 4.95792131

Biểu đồ tỷ lệ đột quỵ theo cao huyết áp

stroke_rate <- aggregate(
  stroke ~ hypertension,
  data = df,
  FUN = mean
)

stroke_rate$percent <- stroke_rate$stroke * 100

bp <- barplot(
  stroke_rate$percent,
  names.arg = c(
    "Không cao huyết áp",
    "Có cao huyết áp"
  ),
  main = "Tỷ lệ đột quỵ theo tình trạng cao huyết áp",
  xlab = "Tình trạng cao huyết áp",
  ylab = "Tỷ lệ đột quỵ (%)",
  col = c("deepskyblue", "tomato"),
  border = "gray30",
  ylim = c(
    0,
    max(stroke_rate$percent) * 1.3
  )
)

text(
  bp,
  stroke_rate$percent,
  labels = paste0(
    round(stroke_rate$percent, 1),
    "%"
  ),
  pos = 3
)

1.4 Mối liên quan giữa tuổi và đột quỵ

modelLog_Age <- glm(
  stroke ~ age,
  family = binomial,
  data = df
)

summary(modelLog_Age)

## 
## Call:
## glm(formula = stroke ~ age, family = binomial, data = df)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -7.231441   0.334980  -21.59   <2e-16 ***
## age          0.074727   0.004922   15.18   <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: 1990.4  on 5109  degrees of freedom
## Residual deviance: 1616.3  on 5108  degrees of freedom
## AIC: 1620.3
## 
## Number of Fisher Scoring iterations: 7

coef(modelLog_Age)

## (Intercept)         age 
## -7.23144139  0.07472671

exp(coef(modelLog_Age))

##  (Intercept)          age 
## 0.0007234773 1.0775896112

exp(confint(modelLog_Age))

## Waiting for profiling to be done...

##                    2.5 %      97.5 %
## (Intercept) 0.0003656084 0.001360584
## age         1.0675366655 1.088350041

Biểu đồ xác suất đột quỵ theo tuổi

age_new <- data.frame(
  age = seq(
    min(df$age, na.rm = TRUE),
    max(df$age, na.rm = TRUE),
    length.out = 200
  )
)

age_new$prob <- predict(
  modelLog_Age,
  newdata = age_new,
  type = "response"
)

plot(
  jitter(df$age),
  jitter(df$stroke, amount = 0.03),
  pch = 19,
  col = "gray60",
  xlab = "Tuổi",
  ylab = "Xác suất đột quỵ",
  main = "Mối liên quan giữa tuổi và đột quỵ",
  ylim = c(0, 1)
)

lines(
  age_new$age,
  age_new$prob,
  col = "red",
  lwd = 3
)

1.5 Mối liên quan giữa bệnh tim và đột quỵ

modelLog_HeartDisease <- glm(
  stroke ~ heart.disease,
  family = binomial,
  data = df
)

summary(modelLog_HeartDisease)

## 
## Call:
## glm(formula = stroke ~ heart.disease, family = binomial, data = df)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -3.13248    0.07188 -43.581   <2e-16 ***
## heart.disease  1.54890    0.17553   8.824   <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: 1990.4  on 5109  degrees of freedom
## Residual deviance: 1930.1  on 5108  degrees of freedom
## AIC: 1934.1
## 
## Number of Fisher Scoring iterations: 6

coef(modelLog_HeartDisease)

##   (Intercept) heart.disease 
##     -3.132476      1.548902

exp(coef(modelLog_HeartDisease))

##   (Intercept) heart.disease 
##    0.04360967    4.70629945

exp(confint(modelLog_HeartDisease))

## Waiting for profiling to be done...

##                   2.5 %     97.5 %
## (Intercept)   0.0377615 0.05005859
## heart.disease 3.3055653 6.58706368

Biểu đồ tỷ lệ đột quỵ theo bệnh tim

stroke_rate_hd <- aggregate(
  stroke ~ heart.disease,
  data = df,
  FUN = mean
)

stroke_rate_hd$percent <- stroke_rate_hd$stroke * 100

bp_hd <- barplot(
  stroke_rate_hd$percent,
  names.arg = c(
    "Không bệnh tim",
    "Có bệnh tim"
  ),
  main = "Tỷ lệ đột quỵ theo tình trạng bệnh tim",
  xlab = "Tình trạng bệnh tim",
  ylab = "Tỷ lệ đột quỵ (%)",
  col = c("deepskyblue", "tomato"),
  border = "gray30",
  ylim = c(
    0,
    max(stroke_rate_hd$percent) * 1.3
  )
)

text(
  bp_hd,
  stroke_rate_hd$percent,
  labels = paste0(
    round(stroke_rate_hd$percent, 1),
    "%"
  ),
  pos = 3
)

1.6 Mô hình logistic đa biến

1.6.1 Hiệu chỉnh theo tuổi, cao huyết áp và bệnh tim

modelLog_adjust <- glm(
  stroke ~ hypertension + age + heart.disease,
  family = binomial,
  data = df
)

summary(modelLog_adjust)

## 
## Call:
## glm(formula = stroke ~ hypertension + age + heart.disease, family = binomial, 
##     data = df)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -7.103034   0.336442 -21.112  < 2e-16 ***
## hypertension   0.457703   0.160691   2.848  0.00439 ** 
## age            0.070354   0.005082  13.844  < 2e-16 ***
## heart.disease  0.408278   0.185995   2.195  0.02816 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1990.4  on 5109  degrees of freedom
## Residual deviance: 1603.6  on 5106  degrees of freedom
## AIC: 1611.6
## 
## Number of Fisher Scoring iterations: 7

coef(modelLog_adjust)

##   (Intercept)  hypertension           age heart.disease 
##   -7.10303358    0.45770288    0.07035408    0.40827777

exp(coef(modelLog_adjust))

##   (Intercept)  hypertension           age heart.disease 
##  0.0008226057  1.5804393552  1.0728879989  1.5042249318

exp(confint(modelLog_adjust))

## Waiting for profiling to be done...

##                      2.5 %      97.5 %
## (Intercept)   0.0004144428 0.001551165
## hypertension  1.1470459811 2.155195462
## age           1.0625296031 1.083924328
## heart.disease 1.0358106250 2.150149029

Biểu đồ Odds Ratio hiệu chỉnh

OR <- exp(coef(modelLog_adjust))[-1]

CI <- exp(confint(modelLog_adjust))[-1, , drop = FALSE]

## Waiting for profiling to be done...

labels <- c(
  "Cao huyết áp",
  "Tuổi",
  "Bệnh tim"
)

y_pos <- 1:3

par(
  mar = c(5, 10, 4, 5)
)

plot(
  OR,
  y_pos,
  xlim = c(
    0.9,
    max(CI[, 2]) * 1.25
  ),
  ylim = c(0.5, 3.5),
  pch = 19,
  cex = 1.5,
  yaxt = "n",
  xlab = "Odds Ratio hiệu chỉnh (OR)",
  ylab = "",
  main = "OR hiệu chỉnh của các yếu tố liên quan đến đột quỵ"
)

axis(
  2,
  at = y_pos,
  labels = labels,
  las = 1
)

segments(
  CI[, 1],
  y_pos,
  CI[, 2],
  y_pos,
  lwd = 2
)

abline(
  v = 1,
  col = "red",
  lty = 2,
  lwd = 2
)

text(
  max(CI[, 2]) * 1.05,
  y_pos,
  labels = paste0(
    "OR = ",
    round(OR, 2),
    " (",
    round(CI[, 1], 2),
    "–",
    round(CI[, 2], 2),
    ")"
  ),
  adj = 0,
  cex = 0.9
)

1.6.2 Thay đổi OR của cao huyết áp sau hiệu chỉnh

OR_Hypertension <- exp(
  coef(modelLog_Hypertension)["hypertension"]
)

OR_Hypertension_adjusted <- exp(
  coef(modelLog_adjust)["hypertension"]
)

Change_Hypertension <- (
  OR_Hypertension_adjusted - OR_Hypertension
) / OR_Hypertension * 100

OR_Hypertension

## hypertension 
##     3.697556

OR_Hypertension_adjusted

## hypertension 
##     1.580439

Change_Hypertension

## hypertension 
##    -57.25719

1.6.3 Thay đổi OR của bệnh tim sau hiệu chỉnh

OR_Heartdisease <- exp(
  coef(modelLog_HeartDisease)["heart.disease"]
)

OR_Heartdisease_adjusted <- exp(
  coef(modelLog_adjust)["heart.disease"]
)

Change_Heartdisease <- (
  OR_Heartdisease_adjusted - OR_Heartdisease
) / OR_Heartdisease * 100

OR_Heartdisease

## heart.disease 
##      4.706299

OR_Heartdisease_adjusted

## heart.disease 
##      1.504225

Change_Heartdisease

## heart.disease 
##     -68.03805

2. Hồi qui logistic – lựa chọn mô hình tối ưu

Chuẩn bị các biến dự báo tiềm năng

predictor <- df[, c(
  "gender",
  "age",
  "hypertension",
  "heart.disease",
  "ever.married",
  "work.type",
  "Residence.type",
  "glucose.level",
  "bmi",
  "smoking"
)]

head(predictor)

##   gender age hypertension heart.disease ever.married work.type Residence.type
## 1 Female  17            0             0           No   Private          Urban
## 2 Female  13            0             0           No  children          Rural
## 3   Male  55            0             0          Yes   Private          Urban
## 4 Female  42            0             0           No   Private          Urban
## 5 Female  31            0             0           No   Private          Urban
## 6 Female  38            0             0          Yes   Private          Urban
##   glucose.level  bmi         smoking
## 1         92.97   NA formerly smoked
## 2         85.81 18.6         Unknown
## 3         89.17 31.5    never smoked
## 4         98.53 18.5    never smoked
## 5        108.89 52.3         Unknown
## 6         91.44   NA         Unknown

dim(predictor)

## [1] 5110   10

str(predictor)

## 'data.frame':    5110 obs. of  10 variables:
##  $ gender        : chr  "Female" "Female" "Male" "Female" ...
##  $ age           : num  17 13 55 42 31 38 24 80 33 20 ...
##  $ hypertension  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ heart.disease : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ ever.married  : chr  "No" "No" "Yes" "No" ...
##  $ work.type     : chr  "Private" "children" "Private" "Private" ...
##  $ Residence.type: chr  "Urban" "Rural" "Urban" "Urban" ...
##  $ glucose.level : num  93 85.8 89.2 98.5 108.9 ...
##  $ bmi           : num  NA 18.6 31.5 18.5 52.3 NA 26.2 NA 42.2 28.8 ...
##  $ smoking       : chr  "formerly smoked" "Unknown" "never smoked" "never smoked" ...

2.1 Bayesian Model Averaging – BMA

library(BMA)

## Loading required package: survival

## 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)

data_bma <- na.omit(
  cbind(
    stroke = df$stroke,
    predictor
  )
)

modelLog_bma <- bic.glm(
  stroke ~ .,
  data = data_bma,
  glm.family = binomial
)

summary(modelLog_bma)

## 
## Call:
## bic.glm.formula(f = stroke ~ ., data = data_bma, glm.family = binomial)
## 
## 
##   10  models were selected
##  Best  5  models (cumulative posterior probability =  0.8293 ): 
## 
##                         p!=0    EV        SD        model 1     model 2   
## Intercept               100    -7.853912  0.406312  -7.767e+00  -7.864e+00
## genderMale                0.0   0.000000  0.000000       .           .    
## age                     100.0   0.071008  0.005774   6.957e-02   7.184e-02
## hypertension             62.8   0.345021  0.298619   5.470e-01       .    
## heart.disease             6.2   0.025457  0.111169       .           .    
## ever.marriedYes           0.0   0.000000  0.000000       .           .    
## work.typeGovt_job         0.0   0.000000  0.000000       .           .    
## work.typeNever_worked     0.0   0.000000  0.000000       .           .    
## work.typePrivate          9.4   0.031162  0.107661       .           .    
## work.typeSelf-employed   12.5  -0.049802  0.145731       .           .    
## Residence.typeUrban       0.0   0.000000  0.000000       .           .    
## glucose.level           100.0   0.005226  0.001270   5.047e-03   5.597e-03
## bmi                       0.0   0.000000  0.000000       .           .    
## smokingnever smoked       0.0   0.000000  0.000000       .           .    
## smokingsmokes             9.2   0.039477  0.137336       .           .    
## smokingUnknown            0.0   0.000000  0.000000       .           .    
##                                                                           
## nVar                                                   3           2      
## BIC                                                 -4.031e+04  -4.031e+04
## post prob                                            0.390       0.237    
##                         model 3     model 4     model 5   
## Intercept               -7.870e+00  -8.079e+00  -7.942e+00
## genderMale                   .           .           .    
## age                      7.290e-02   7.143e-02   7.107e-02
## hypertension             5.617e-01   5.533e-01   5.463e-01
## heart.disease                .           .           .    
## ever.marriedYes              .           .           .    
## work.typeGovt_job            .           .           .    
## work.typeNever_worked        .           .           .    
## work.typePrivate             .       3.343e-01       .    
## work.typeSelf-employed  -4.040e-01       .           .    
## Residence.typeUrban          .           .           .    
## glucose.level            4.974e-03   5.003e-03   5.088e-03
## bmi                          .           .           .    
## smokingnever smoked          .           .           .    
## smokingsmokes                .           .       4.308e-01
## smokingUnknown               .           .           .    
##                                                           
## nVar                       4           4           4      
## BIC                     -4.031e+04  -4.030e+04  -4.030e+04
## post prob                0.085       0.061       0.057    
## 
##   1  observations deleted due to missingness.

Dự báo bằng BMA

pred_bma <- predict(
  modelLog_bma,
  newdata = data_bma
)

head(pred_bma)

##           2           3           4           5           7           9 
## 0.001534093 0.030783386 0.013069971 0.006359425 0.003655871 0.006535131

class_bma <- ifelse(
  pred_bma >= 0.5,
  1,
  0
)

table(
  Actual = data_bma$stroke,
  Predicted = class_bma
)

##       Predicted
## Actual    0
##      0 4700
##      1  209

accuracy_bma <- mean(
  class_bma == data_bma$stroke
)

accuracy_bma

## [1] 0.9574251

2.2 LASSO

library(glmnet)

## Loading required package: Matrix

## Loaded glmnet 5.0

library(Matrix)

data_lasso <- na.omit(
  cbind(
    stroke = df$stroke,
    predictor
  )
)

Y <- data_lasso$stroke

X <- model.matrix(
  stroke ~ .,
  data = data_lasso
)[, -1]

set.seed(123)

modelLog_lasso <- cv.glmnet(
  x = X,
  y = Y,
  family = "binomial",
  alpha = 1
)

modelLog_lasso$lambda.min

## [1] 0.001500674

coef(
  modelLog_lasso,
  s = "lambda.min"
)

## 16 x 1 sparse Matrix of class "dgCMatrix"
##                          lambda.min
## (Intercept)            -7.522265374
## genderMale              .          
## age                     0.066251995
## hypertension            0.504050399
## heart.disease           0.351621921
## ever.marriedYes         .          
## work.typeGovt_job       .          
## work.typeNever_worked   .          
## work.typePrivate        0.084877598
## work.typeSelf-employed -0.199300329
## Residence.typeUrban     .          
## glucose.level           0.004446046
## bmi                     .          
## smokingnever smoked     .          
## smokingsmokes           0.245723178
## smokingUnknown         -0.110511985

Dự báo bằng LASSO

pred_lasso <- predict(
  modelLog_lasso,
  newx = X,
  s = "lambda.min",
  type = "response"
)

pred_lasso <- as.vector(pred_lasso)

head(pred_lasso)

## [1] 0.001675467 0.032387031 0.014532408 0.006626974 0.004435683 0.007657374

class_lasso <- ifelse(
  pred_lasso >= 0.5,
  1,
  0
)

table(
  Actual = Y,
  Predicted = class_lasso
)

##       Predicted
## Actual    0    1
##      0 4700    0
##      1  208    1

accuracy_lasso <- mean(
  class_lasso == Y
)

accuracy_lasso

## [1] 0.9576288

So sánh BMA và LASSO bằng ROC-AUC

library(pROC)

## Type 'citation("pROC")' for a citation.

## 
## Attaching package: 'pROC'

## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var

roc_bma <- roc(
  data_bma$stroke,
  pred_bma
)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

roc_lasso <- roc(
  Y,
  pred_lasso
)

## Setting levels: control = 0, case = 1
## Setting direction: controls < cases

par(
  mfrow = c(2, 2),
  mar = c(4, 4, 3, 1)
)

imageplot.bma(modelLog_bma)

plot(
  modelLog_lasso$glmnet.fit,
  xvar = "lambda"
)

plot(
  roc_bma,
  main = paste(
    "BMA ROC - AUC =",
    round(auc(roc_bma), 3)
  )
)

plot(
  roc_lasso,
  main = paste(
    "LASSO ROC - AUC =",
    round(auc(roc_lasso), 3)
  )
)

par(mfrow = c(1, 1))

3. Hồi qui logistic – đánh giá mô hình

dbf <- df[, c(
  "stroke",
  "age",
  "hypertension",
  "heart.disease",
  "glucose.level"
)]

dbf <- na.omit(dbf)

summary(dbf)

##      stroke             age         hypertension     heart.disease    
##  Min.   :0.00000   Min.   : 0.08   Min.   :0.00000   Min.   :0.00000  
##  1st Qu.:0.00000   1st Qu.:25.00   1st Qu.:0.00000   1st Qu.:0.00000  
##  Median :0.00000   Median :45.00   Median :0.00000   Median :0.00000  
##  Mean   :0.04873   Mean   :43.23   Mean   :0.09746   Mean   :0.05401  
##  3rd Qu.:0.00000   3rd Qu.:61.00   3rd Qu.:0.00000   3rd Qu.:0.00000  
##  Max.   :1.00000   Max.   :82.00   Max.   :1.00000   Max.   :1.00000  
##  glucose.level   
##  Min.   : 55.12  
##  1st Qu.: 77.25  
##  Median : 91.89  
##  Mean   :106.15  
##  3rd Qu.:114.09  
##  Max.   :271.74

3.1 Đánh giá khả năng phân định bằng ROC và AUC

model_Log_dbf <- glm(
  stroke ~ age + hypertension + heart.disease + glucose.level,
  family = binomial,
  data = dbf
)

summary(model_Log_dbf)

## 
## Call:
## glm(formula = stroke ~ age + hypertension + heart.disease + glucose.level, 
##     family = binomial, data = dbf)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -7.489396   0.357879 -20.927  < 2e-16 ***
## age            0.068926   0.005140  13.410  < 2e-16 ***
## hypertension   0.381410   0.162599   2.346  0.01899 *  
## heart.disease  0.329965   0.187724   1.758  0.07880 .  
## glucose.level  0.004121   0.001162   3.546  0.00039 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1990.4  on 5109  degrees of freedom
## Residual deviance: 1591.5  on 5105  degrees of freedom
## AIC: 1601.5
## 
## Number of Fisher Scoring iterations: 7

predict_Log_dbf <- predict(
  model_Log_dbf,
  type = "response"
)

head(predict_Log_dbf)

##           1           2           3           4           5           6 
## 0.002639410 0.001946519 0.034521770 0.014942364 0.007362297 0.011058686

roc_Log_dbf <- roc(
  dbf$stroke,
  predict_Log_dbf
)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

plot(
  roc_Log_dbf,
  main = paste(
    "ROC Curve - AUC =",
    round(auc(roc_Log_dbf), 3)
  )
)

auc(roc_Log_dbf)

## Area under the curve: 0.8442

Nhận xét

Khả năng phân định của mô hình được đánh giá thông qua đường cong ROC và diện tích dưới đường cong AUC. AUC càng gần 1 cho thấy mô hình càng có khả năng phân biệt tốt giữa nhóm có đột quỵ và nhóm không đột quỵ.

Nếu AUC khoảng 0.8 trở lên, mô hình có khả năng phân định tốt.

3.2 Đánh giá mô hình bằng Bootstrap

library(rms)

## Loading required package: Hmisc

## 
## Attaching package: 'Hmisc'

## The following objects are masked from 'package:table1':
## 
##     label, label<-, units

## The following objects are masked from 'package:base':
## 
##     format.pval, units

model_rms <- lrm(
  stroke ~ age + hypertension + heart.disease + glucose.level,
  data = dbf,
  x = TRUE,
  y = TRUE
)

set.seed(123)

boot_dbf <- validate(
  model_rms,
  method = "boot",
  B = 100
)

boot_dbf

##           index.orig training    test optimism index.corrected   Lower  Upper
## Dxy           0.6884   0.6893  0.6868   0.0025          0.6859  0.6436 0.7252
## R2            0.2328   0.2339  0.2310   0.0028          0.2299  0.1936 0.2610
## Intercept     0.0000   0.0000 -0.0117   0.0117         -0.0117 -0.3414 0.3426
## Slope         1.0000   1.0000  0.9895   0.0105          0.9895  0.8465 1.1271
## Emax          0.0000   0.0000  0.0322  -0.0322          0.0322 -0.0116 0.1064
## D             0.0779   0.0779  0.0773   0.0007          0.0772  0.0632 0.0906
## U            -0.0004  -0.0004  0.0000  -0.0004          0.0000 -0.0006 0.0012
## Q             0.0783   0.0783  0.0773   0.0010          0.0772  0.0629 0.0905
## B             0.0421   0.0418  0.0422  -0.0004          0.0425  0.0378 0.0471
## g             1.9252   1.9330  1.9090   0.0240          1.9012  1.6578 2.1038
## gp            0.0633   0.0630  0.0630  -0.0001          0.0633  0.0556 0.0709
##             n
## Dxy       100
## R2        100
## Intercept 100
## Slope     100
## Emax      100
## D         100
## U         100
## Q         100
## B         100
## g         100
## gp        100

AUC hiệu chỉnh bằng Bootstrap

Dxy_boot <- boot_dbf[
  "Dxy",
  "index.corrected"
]

AUC_boot <- (
  Dxy_boot + 1
) / 2

AUC_boot

## [1] 0.8429432

Nhận xét

Bootstrap được sử dụng để đánh giá độ ổn định của mô hình và hiệu chỉnh mức độ lạc quan (optimism) khi mô hình được đánh giá trên chính dữ liệu dùng để xây dựng.

Có thể so sánh:

auc(roc_Log_dbf)

## Area under the curve: 0.8442

AUC_boot

## [1] 0.8429432

Nếu AUC sau bootstrap gần với AUC ban đầu, mô hình tương đối ổn định và ít có dấu hiệu overfitting.

3.3 Đánh giá mô hình bằng phương pháp chia dữ liệu với caret

library(caret)

## Loading required package: ggplot2

## Loading required package: lattice

## Registered S3 method overwritten by 'plyr':
##   method    from  
##   [.indexed table1

## 
## Attaching package: 'caret'

## The following object is masked from 'package:survival':
## 
##     cluster

caret_data <- dbf

caret_data$stroke <- factor(
  caret_data$stroke,
  levels = c(0, 1),
  labels = c("No", "Yes")
)

set.seed(123)

index <- createDataPartition(
  caret_data$stroke,
  p = 0.7,
  list = FALSE
)

train_data <- caret_data[index, ]

test_data <- caret_data[-index, ]

Xây dựng mô hình trên tập train

model_caret <- train(
  stroke ~ age + hypertension + heart.disease + glucose.level,
  data = train_data,
  method = "glm",
  family = binomial
)

model_caret

## Generalized Linear Model 
## 
## 3578 samples
##    4 predictor
##    2 classes: 'No', 'Yes' 
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 3578, 3578, 3578, 3578, 3578, 3578, ... 
## Resampling results:
## 
##   Accuracy   Kappa      
##   0.9528566  0.001178986

Dự báo trên tập test

class_caret <- predict(
  model_caret,
  newdata = test_data
)

confusionMatrix(
  class_caret,
  test_data$stroke,
  positive = "Yes"
)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   No  Yes
##        No  1458   74
##        Yes    0    0
##                                           
##                Accuracy : 0.9517          
##                  95% CI : (0.9397, 0.9619)
##     No Information Rate : 0.9517          
##     P-Value [Acc > NIR] : 0.5309          
##                                           
##                   Kappa : 0               
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.0000          
##             Specificity : 1.0000          
##          Pos Pred Value :    NaN          
##          Neg Pred Value : 0.9517          
##              Prevalence : 0.0483          
##          Detection Rate : 0.0000          
##    Detection Prevalence : 0.0000          
##       Balanced Accuracy : 0.5000          
##                                           
##        'Positive' Class : Yes             
##

ROC và AUC trên tập test

prob_caret <- predict(
  model_caret,
  newdata = test_data,
  type = "prob"
)[, "Yes"]

roc_caret <- roc(
  test_data$stroke,
  prob_caret
)

## Setting levels: control = No, case = Yes

## Setting direction: controls < cases

auc_caret <- auc(roc_caret)

auc_caret

## Area under the curve: 0.8179

plot(
  roc_caret,
  main = paste(
    "ROC Test Set - AUC =",
    round(auc_caret, 3)
  )
)

Nhận xét

Phương pháp chia dữ liệu giúp đánh giá khả năng dự báo của mô hình trên các quan sát không được sử dụng trong quá trình xây dựng mô hình.

Các chỉ số trong confusionMatrix() như Accuracy, Sensitivity và Specificity cùng với AUC trên tập test được sử dụng để đánh giá hiệu quả dự báo.

AUC trên tập test càng cao và càng gần với AUC của mô hình ban đầu thì mô hình càng có khả năng tổng quát hóa tốt.

Kết luận

Phân tích hồi quy logistic cho phép đánh giá mối liên quan giữa các yếu tố nguy cơ và khả năng xảy ra đột quỵ. Các phương pháp BMA và LASSO được sử dụng để hỗ trợ lựa chọn mô hình và biến dự báo.

Khả năng dự báo của mô hình được đánh giá thông qua ba hướng chính:

ROC và AUC để đánh giá khả năng phân định.
Bootstrap để đánh giá độ ổn định và hiệu chỉnh optimism.
Chia dữ liệu train-test để đánh giá khả năng dự báo trên dữ liệu độc lập.

Việc kết hợp các phương pháp này giúp đánh giá mô hình logistic toàn diện hơn về cả khả năng phân định, độ ổn định và khả năng tổng quát hóa.

Hồi quy Logistic trong dự báo và đánh giá nguy cơ đột quỵ

Chính Dương

2026-08-18

1. Hồi qui logistic – hiệu chỉnh yếu tố nhiễu

1.1 Đọc dữ liệu

1.2 Mô tả dữ liệu theo tình trạng đột quỵ

1.3 Mối liên quan giữa cao huyết áp và đột quỵ

Biểu đồ tỷ lệ đột quỵ theo cao huyết áp

1.4 Mối liên quan giữa tuổi và đột quỵ

Biểu đồ xác suất đột quỵ theo tuổi

1.5 Mối liên quan giữa bệnh tim và đột quỵ

Biểu đồ tỷ lệ đột quỵ theo bệnh tim

1.6 Mô hình logistic đa biến

1.6.1 Hiệu chỉnh theo tuổi, cao huyết áp và bệnh tim

Biểu đồ Odds Ratio hiệu chỉnh

1.6.2 Thay đổi OR của cao huyết áp sau hiệu chỉnh

1.6.3 Thay đổi OR của bệnh tim sau hiệu chỉnh

2. Hồi qui logistic – lựa chọn mô hình tối ưu

Chuẩn bị các biến dự báo tiềm năng

2.1 Bayesian Model Averaging – BMA

Dự báo bằng BMA

2.2 LASSO

Dự báo bằng LASSO

So sánh BMA và LASSO bằng ROC-AUC

3. Hồi qui logistic – đánh giá mô hình

3.1 Đánh giá khả năng phân định bằng ROC và AUC

Nhận xét

3.2 Đánh giá mô hình bằng Bootstrap

AUC hiệu chỉnh bằng Bootstrap

Nhận xét

3.3 Đánh giá mô hình bằng phương pháp chia dữ liệu với caret

Xây dựng mô hình trên tập train

Dự báo trên tập test

ROC và AUC trên tập test

Nhận xét

Kết luận