D5_exercise
Ha Tan Duc
2022-09-18
Việc 1:
Đọc dữ liệu nghiên cứu tim mạch Framingham dataset.csv
và gọi đối tượng fmh. Có những biến số liên quan đến bệnh
tim mạch. Các biến liên quan:
age tuổi \ sex giới tính (1/2, 1=male,
2=female) \ sysbp systolic blood pressure \
diasbp diastolic blood pressure \ tot.chol
Total cholesterol \ heart.rate nhịp tim \ bmi
tỉ trọng cơ thể \ smoker hút thuốc lá (0/1) \
diabetes tiểu đường (0/1) \ stroke đột quị
(0/1) \ cvd cardiovascular disease (0/1) \
death tử vong (0/1) \
fmh = read.csv("D:\\OneDrive\\Statistical courses\\Can Tho University of Medicine_Sep2022\\Data for practice\\Framingham dataset.csv")
names (fmh)## [1] "id" "sex" "tot.chol" "age" "sysbp"
## [6] "diasbp" "smoker" "cigs.day" "bmi" "diabetes"
## [11] "bpmed" "heart.rate" "glucose" "educ" "prev.chd"
## [16] "prev.ap" "prev.mi" "prev.stroke" "prev.hyp" "time"
## [21] "period" "hdlc" "ldlc" "death" "angina"
## [26] "hosp.mi" "mi.fchd" "any.chd" "stroke" "cvd"
## [31] "hypertension" "time.ap" "time.mi" "time.mi.1" "time.chd"
## [36] "time.stroke" "time.cvd" "time.dth" "time.hyp"head (fmh)## id sex tot.chol age sysbp diasbp smoker cigs.day bmi diabetes bpmed
## 1 2448 1 195 39 106.0 70.0 0 0 26.97 0 0
## 2 2448 1 209 52 121.0 66.0 0 0 NA 0 0
## 3 6238 2 250 46 121.0 81.0 0 0 28.73 0 0
## 4 6238 2 260 52 105.0 69.5 0 0 29.43 0 0
## 5 6238 2 237 58 108.0 66.0 0 0 28.50 0 0
## 6 9428 1 245 48 127.5 80.0 1 20 25.34 0 0
## heart.rate glucose educ prev.chd prev.ap prev.mi prev.stroke prev.hyp time
## 1 80 77 4 0 0 0 0 0 0
## 2 69 92 4 0 0 0 0 0 4628
## 3 95 76 2 0 0 0 0 0 0
## 4 80 86 2 0 0 0 0 0 2156
## 5 80 71 2 0 0 0 0 0 4344
## 6 75 70 1 0 0 0 0 0 0
## period hdlc ldlc death angina hosp.mi mi.fchd any.chd stroke cvd hypertension
## 1 1 NA NA 0 0 1 1 1 0 1 0
## 2 3 31 178 0 0 1 1 1 0 1 0
## 3 1 NA NA 0 0 0 0 0 0 0 0
## 4 2 NA NA 0 0 0 0 0 0 0 0
## 5 3 54 141 0 0 0 0 0 0 0 0
## 6 1 NA NA 0 0 0 0 0 0 0 0
## time.ap time.mi time.mi.1 time.chd time.stroke time.cvd time.dth time.hyp
## 1 8766 6438 6438 6438 8766 6438 8766 8766
## 2 8766 6438 6438 6438 8766 6438 8766 8766
## 3 8766 8766 8766 8766 8766 8766 8766 8766
## 4 8766 8766 8766 8766 8766 8766 8766 8766
## 5 8766 8766 8766 8766 8766 8766 8766 8766
## 6 8766 8766 8766 8766 8766 8766 8766 8766Tạo ra biến mới gender từ biến sex.
fmh$gender [fmh$sex == 1] = "Male"
fmh$gender [fmh$sex == 2] = "Female"Chúng ta sẽ tập trung vào phân tích biến số stroke (đột
quị). Các biến tiên lượng hay “yếu tố nguy cơ” bao gồm age,
gender, sysbp, tot.chol,
bmi, smoker, diabetes.
table1::table1 (~ age + gender + sysbp + tot.chol + bmi + factor(smoker) + factor(diabetes)| stroke, data = fmh)## Warning in table1.formula(~age + gender + sysbp + tot.chol + bmi +
## factor(smoker) + : Terms to the right of '|' in formula 'x' define table columns
## and are expected to be factors with meaningful labels.| 0 (N=10566) |
1 (N=1061) |
Overall (N=11627) |
|
|---|---|---|---|
| age | |||
| Mean (SD) | 54.2 (9.45) | 60.4 (8.90) | 54.8 (9.56) |
| Median [Min, Max] | 54.0 [32.0, 81.0] | 61.0 [34.0, 81.0] | 54.0 [32.0, 81.0] |
| gender | |||
| Female | 6032 (57.1%) | 573 (54.0%) | 6605 (56.8%) |
| Male | 4534 (42.9%) | 488 (46.0%) | 5022 (43.2%) |
| sysbp | |||
| Mean (SD) | 135 (21.8) | 151 (26.9) | 136 (22.8) |
| Median [Min, Max] | 131 [83.5, 244] | 147 [94.0, 295] | 132 [83.5, 295] |
| tot.chol | |||
| Mean (SD) | 241 (44.9) | 243 (49.8) | 241 (45.4) |
| Median [Min, Max] | 238 [113, 638] | 240 [107, 696] | 238 [107, 696] |
| Missing | 376 (3.6%) | 33 (3.1%) | 409 (3.5%) |
| bmi | |||
| Mean (SD) | 25.8 (4.03) | 26.8 (4.68) | 25.9 (4.10) |
| Median [Min, Max] | 25.4 [14.4, 56.8] | 26.4 [16.9, 55.3] | 25.5 [14.4, 56.8] |
| Missing | 42 (0.4%) | 10 (0.9%) | 52 (0.4%) |
| factor(smoker) | |||
| 0 | 5981 (56.6%) | 617 (58.2%) | 6598 (56.7%) |
| 1 | 4585 (43.4%) | 444 (41.8%) | 5029 (43.3%) |
| factor(diabetes) | |||
| 0 | 10164 (96.2%) | 933 (87.9%) | 11097 (95.4%) |
| 1 | 402 (3.8%) | 128 (12.1%) | 530 (4.6%) |
Dùng ggplot để thể hiện mối liên quan giữa mỗi yếu tố nguy cơ và stroke.
library (ggplot2)
ggplot(fmh, aes(x= factor(stroke), y= bmi, fill= factor(stroke))) +
geom_boxplot() + geom_jitter(width=0.25, alpha=0.1)## Warning: Removed 52 rows containing non-finite values (stat_boxplot).## Warning: Removed 52 rows containing missing values (geom_point).
Việc 2:
Mô hình 1 với biến độc lập là biến phân loại. Xem stroke là biến phụ thuộc và gender là yếu tố nguy cơ. Thường, các chuyên gia cho rằng nam giới có nguy cơ đột quị cao hơn nữ. Các bạn hãy dùng dữ liệu này để kiểm định giả thuyết trên.
- Trong số 5022 nam, có bao nhiêu người bị đột quị? Trong số 6605 nữ, bao nhiêu người bị đột quị? Tỉ lệ đột quị ở nam giới có cao hơn nữ giới? Dùng phương pháp Ki bình phương để so sánh 2 tỉ lệ.
DescTools::Desc (factor (fmh$stroke) ~ fmh$gender)## ------------------------------------------------------------------------------
## factor(fmh$stroke) ~ fmh$gender
##
## Summary:
## n: 11'627, rows: 2, columns: 2
##
## Pearson's Chi-squared test (cont. adj):
## X-squared = 3.6107, df = 1, p-value = 0.05741
## Fisher's exact test p-value = 0.05512
## McNemar's chi-squared = 3070.6, df = 1, p-value < 2.2e-16
##
## estimate lwr.ci upr.ci'
##
## odds ratio 1.133 0.998 1.286
## rel. risk (col1) 1.057 0.998 1.120
## rel. risk (col2) 0.933 0.871 0.999
##
##
## Contingency Coeff. 0.018
## Cramer's V 0.018
## Kendall Tau-b 0.018
##
##
## fmh$gender Female Male Sum
## factor(fmh$stroke)
##
## 0 freq 6'032 4'534 10'566
## perc 51.9% 39.0% 90.9%
## p.row 57.1% 42.9% .
## p.col 91.3% 90.3% .
##
## 1 freq 573 488 1'061
## perc 4.9% 4.2% 9.1%
## p.row 54.0% 46.0% .
## p.col 8.7% 9.7% .
##
## Sum freq 6'605 5'022 11'627
## perc 56.8% 43.2% 100.0%
## p.row . . .
## p.col . . .
##
##
## ----------
## ' 95% conf. level
- Các bạn hãy dùng mô hình hồi qui logistic để ước tính odds đột quị ở nam giới so với nữ giới. Diễn giải kết quả phân tích. Sự khác biệt có ý nghĩa thống kê không?
summary (glm (stroke ~ gender, family = binomial, data= fmh))##
## Call:
## glm(formula = stroke ~ gender, family = binomial, data = fmh)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.4522 -0.4522 -0.4260 -0.4260 2.2112
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.35395 0.04371 -53.848 <2e-16 ***
## genderMale 0.12490 0.06466 1.932 0.0534 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 7102.4 on 11626 degrees of freedom
## Residual deviance: 7098.7 on 11625 degrees of freedom
## AIC: 7102.7
##
## Number of Fisher Scoring iterations: 5fmh$gender = as.factor (fmh$gender)
m= glm (stroke ~ gender, family = binomial, data= fmh)
epiDisplay::logistic.display(m)##
## Logistic regression predicting stroke
##
## OR(95%CI) P(Wald's test) P(LR-test)
## gender: Male vs Female 1.13 (1,1.29) 0.053 0.054
##
## Log-likelihood = -3549.346
## No. of observations = 11627
## AIC value = 7102.6921Việc 3:
Mô hình 2 với biến độc lập là biến liên tục. Y văn chỉ ra rằng người
bị đột quị thường có huyết áp cao. Hãy kiểm tra mối liên quan này qua dữ
liệu Framingham (stroke là biến phụ thuộc, và
sysbp là yếu tố nguy cơ).
- Vẽ biểu đồ hộp thể hiện mối liên quan giữa
strokevàsysbp. Các bạn có thấysysbpở nhóm đột quị cao hơn nhóm không bị đột quị?
library (ggplot2)
ggplot(fmh, aes(x= factor(stroke), y= sysbp, fill= factor(stroke))) +
geom_boxplot() + geom_jitter(width=0.25, alpha=0.1)
- Sử dụng kiểm định t để kiểm tra giả thuyết rằng
sysbpcao hơn ở người bị đột quị so với người không bị đột quị, báo cáo khoảng tin cậy 95% và trị số P.
hist (fmh$sysbp);
t.test(fmh$sysbp ~ factor(fmh$stroke))##
## Welch Two Sample t-test
##
## data: fmh$sysbp by factor(fmh$stroke)
## t = -18.931, df = 1203.6, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -17.83470 -14.48521
## sample estimates:
## mean in group 0 mean in group 1
## 134.8495 151.0094- Dùng mô hình hồi qui logistic để ước tính odds đột quị cho mỗi 10 mmHg tăng huyết áp. Viết một dòng mô tả kết quả và diễn giải kết quả.
fmh$sys = fmh$sysbp/10
summary (glm (stroke ~ sys, family = binomial, data= fmh));##
## Call:
## glm(formula = stroke ~ sys, family = binomial, data = fmh)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3440 -0.4530 -0.3701 -0.3094 2.6820
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.04551 0.19005 -31.81 <2e-16 ***
## sys 0.26348 0.01265 20.83 <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: 7102.4 on 11626 degrees of freedom
## Residual deviance: 6676.4 on 11625 degrees of freedom
## AIC: 6680.4
##
## Number of Fisher Scoring iterations: 5epiDisplay::logistic.display(glm (stroke ~ sys, family = binomial, data= fmh))##
## Logistic regression predicting stroke
##
## OR(95%CI) P(Wald's test) P(LR-test)
## sys (cont. var.) 1.3 (1.27,1.33) < 0.001 < 0.001
##
## Log-likelihood = -3338.1754
## No. of observations = 11627
## AIC value = 6680.3508Việc 4: Mô hình 2 biến và hiệu chỉnh
- Đánh giá ảnh hưởng của hút thuốc lá (smoker) theo mô hình sau đây:
stroke ~ α + β*smokerCác bạn chú ý rằng kết quả cho thấy người hút thuốc lá có nguy cơ đột quị thấp hơn người không hút thuốc lá.
epiDisplay::logistic.display(glm (stroke ~ smoker, family = binomial, data= fmh))##
## Logistic regression predicting stroke
##
## OR(95%CI) P(Wald's test) P(LR-test)
## smoker: 1 vs 0 0.94 (0.83,1.07) 0.332 0.332
##
## Log-likelihood = -3550.735
## No. of observations = 11627
## AIC value = 7105.4699- Đánh giá ảnh hưởng của độ tuổi (age)
epiDisplay::logistic.display(glm (stroke ~ age, family = binomial, data= fmh))##
## Logistic regression predicting stroke
##
## OR(95%CI) P(Wald's test) P(LR-test)
## age (cont. var.) 1.07 (1.06,1.08) < 0.001 < 0.001
##
## Log-likelihood = -3351.0549
## No. of observations = 11627
## AIC value = 6706.1098- So sánh độ tuổi (age) giữa người hút thuốc lá và không hút thuốc lá (qua kiểm định t) và biểu đồ hộp. Có phải người hút thuốc lá có tuổi thấp hơn người không hút thuốc lá?
library (ggplot2)
ggplot(fmh, aes(x= factor(smoker), y= age, fill= factor(smoker))) +
geom_boxplot() + geom_jitter(width=0.25, alpha=0.1)
t.test (fmh$age ~ fmh$smoker)##
## Welch Two Sample t-test
##
## data: fmh$age by fmh$smoker
## t = 28.502, df = 11256, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 4.541912 5.212768
## sample estimates:
## mean in group 0 mean in group 1
## 56.90239 52.02505- Phân tích mô hình sau đây:
stroke ~ α + β*smoker + γ*age
fmh$smoker = as.factor(fmh$smoker)
epiDisplay::logistic.display(glm (stroke ~ smoker + age, family = binomial, data= fmh))##
## Logistic regression predicting stroke
##
## crude OR(95%CI) adj. OR(95%CI) P(Wald's test) P(LR-test)
## smoker: 1 vs 0 0.94 (0.83,1.07) 1.37 (1.2,1.57) < 0.001 < 0.001
##
## age (cont. var.) 1.07 (1.06,1.08) 1.08 (1.07,1.08) < 0.001 < 0.001
##
## Log-likelihood = -3340.7762
## No. of observations = 11627
## AIC value = 6687.5524Việc 5: Tìm mô hình tối ưu
Nghiên cứu Framingham có biến phụ thuộc là stroke. Hãy xem các biến
sau đây là yếu tố nguy cơ: age, gender,
sysbp, tot.chol, bmi,
smoker, diabetes. Hãy dùng phương pháp BMA
(Bayesian Model Averaging) để tìm biến liên quan đến stroke.
library(BMA)## Loading required package: survival## Loading required package: leaps## Warning: package 'leaps' was built under R version 4.0.5## Loading required package: robustbase## Warning: package 'robustbase' was built under R version 4.0.5##
## Attaching package: 'robustbase'## The following object is masked from 'package:survival':
##
## heart## Loading required package: inline## Warning: package 'inline' was built under R version 4.0.5## Loading required package: rrcov## Scalable Robust Estimators with High Breakdown Point (version 1.7-1)fmh$diabetes = as.factor(fmh$diabetes)
xvars = fmh [, c('age', 'gender', 'sysbp', 'tot.chol', 'bmi', 'smoker', 'diabetes')]
r2= bic.glm(x= xvars, y= fmh$stroke, glm.family = binomial, maxCol= 30)## Warning in bic.glm.data.frame(x = xvars, y = fmh$stroke, glm.family =
## binomial, : There were 454 records deleted due to NA'ssummary (r2)##
## Call:
## bic.glm.data.frame(x = xvars, y = fmh$stroke, glm.family = binomial, maxCol = 30)
##
##
## 6 models were selected
## Best 5 models (cumulative posterior probability = 0.9726 ):
##
## p!=0 EV SD model 1 model 2 model 3
## Intercept 100 -8.663e+00 0.4411991 -8.377e+00 -8.970e+00 -8.500e+00
## age 100.0 5.476e-02 0.0041749 5.406e-02 5.560e-02 5.390e-02
## gender 20.9
## .Male 3.838e-02 0.0811847 . . 1.859e-01
## sysbp 100.0 1.913e-02 0.0015413 1.952e-02 1.847e-02 1.995e-02
## tot.chol 5.5 -8.453e-05 0.0003928 . . .
## bmi 47.0 1.145e-02 0.0133367 . 2.443e-02 .
## smoker 100.0
## .1 3.933e-01 0.0761090 3.842e-01 4.180e-01 3.532e-01
## diabetes 100.0
## .1 7.045e-01 0.1190784 7.220e-01 6.914e-01 7.092e-01
##
## nVar 4 5 5
## BIC -9.792e+04 -9.792e+04 -9.792e+04
## post prob 0.381 0.355 0.121
## model 4 model 5
## Intercept -9.081e+00 -8.053e+00
## age 5.543e-02 5.421e-02
## gender
## .Male 1.796e-01 .
## sysbp 1.893e-02 1.980e-02
## tot.chol . -1.522e-03
## bmi 2.395e-02 .
## smoker
## .1 3.871e-01 3.789e-01
## diabetes
## .1 6.800e-01 7.125e-01
##
## nVar 6 5
## BIC -9.791e+04 -9.791e+04
## post prob 0.088 0.028Đánh giá relative importance của các preditors trong mô hình (sử dụng
caret::varImp)
m2= glm (stroke ~ age + smoker + sysbp, family = binomial, data= fmh)
caret::varImp(m2, scale = F)## Overall
## age 14.030811
## smoker1 5.223843
## sysbp 14.682717# Vẽ biểu đồ biểu diễn Odds ratio
library(GGally); library(ggplot2)## Warning: package 'GGally' was built under R version 4.0.5## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2options(scipen= 999) #trình bày số theo quy ước thông thường
ggcoef(m2, exponentiate = TRUE, exclude_intercept=T, mapping = aes(x = estimate, y = term, size = p.value)) + scale_size_continuous(trans = "reverse")