Bài tập ngày 3 + 4
2024-09-10
Phân tích so sánh
Việc 1. Đọc dữ liệu “Obesity data.csv” vào R và gọi dữ liệu là “ob”
ob=read.csv("C:\\Dr Nhung\\Khóa học phân tích dữ liệu trong NCKH bv 108 tháng 9.2024\\obesity data.csv")Việc 2. So sánh biến liên tục cho 2 nhóm
Bài tập ngày 4
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.4.1library(gridExtra) ## Warning: package 'gridExtra' was built under R version 4.4.1p = ggplot(data = ob, aes(x = age))
p1 = p + geom_histogram(fill = "blue", col = "white") + labs(x = "Tuổi (năm)", y = "Số người", title = "Phân bố tuổi")
p = ggplot(data = ob, aes(x = wbbmd))
p2 = p + geom_histogram(fill = "blue", col = "white") + labs(x = "Mật độ xương toàn thân (g/cm2)", y = "Số người", title = "Phân bố MĐX toàn thân")
grid.arrange(p1, p2, ncol = 2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Biểu đồ tán xạ đánh giá mối liên quan giữa tuổi và MĐX toàn thân
p = ggplot(data = ob, aes(x = age, y = wbbmd))
p + geom_point() + geom_smooth(method = "lm", formula = y~ x)
Phân tích tương quan đánh giá mối liên quan giữa tuổi và MĐX toàn thân
Mô hình tuyến tính đánh giá mối liên quan giữa tuổi và MĐX toàn thân
Thực hiện mô hình hồi qui
m.1 = lm(wbbmd ~ age, data = ob)
summary(m.1)##
## Call:
## lm(formula = wbbmd ~ age, data = ob)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32749 -0.07268 -0.00533 0.06793 0.33178
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.1450766 0.0084638 135.29 <2e-16 ***
## age -0.0028914 0.0001686 -17.15 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1015 on 1215 degrees of freedom
## Multiple R-squared: 0.195, Adjusted R-squared: 0.1943
## F-statistic: 294.3 on 1 and 1215 DF, p-value: < 2.2e-16par(mfrow = c(2,2))
plot(m.1)
# ghi chú cho trên: Lệnh par(mfrow = c(2,2)) giúp tạo 4 biểu đồ/1 trang
giấy Residual = phần dư
Phương trình: MDX = 1.145 - 0.003*age
Đánh giá giả định của mô hình
Q-Q plot = Q-Q Residual = đánh giá phân bố (chuẩn hay ko)
Việc 4. Đánh giá mối liên quan giữa tuổi và giới tính với MĐX toàn thân
4.1 Biểu đồ tán xạ đánh giá mối liên quan giữa tuổi với MĐX toàn thân theo giới tính
p = ggplot(data = ob, aes(x = age, y = wbbmd, fill = gender, col = gender))
p1 = p + geom_point() + geom_smooth() + labs(x = "Tuổi (năm)", y = "Mật độ xương toàn thân (g/cm2)") + ggtitle("Liên quan giữa tuổi và MĐX theo giới tính")
p1## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## 4.2 Nhận xét về ảnh hưởng của giới tính lên mối liên quan giữa tuổi
với MĐX toàn thân 4.3 Thực hiện mô hình hồi qui tuyến tính
m.2 = lm(wbbmd ~ age + gender, data = ob)
summary(m.2)##
## Call:
## lm(formula = wbbmd ~ age + gender, data = ob)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.36272 -0.06658 -0.00411 0.06549 0.34473
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.118288 0.008636 129.485 <2e-16 ***
## age -0.002691 0.000164 -16.408 <2e-16 ***
## genderM 0.059417 0.006230 9.537 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09798 on 1214 degrees of freedom
## Multiple R-squared: 0.2511, Adjusted R-squared: 0.2498
## F-statistic: 203.5 on 2 and 1214 DF, p-value: < 2.2e-16Tuổi liên quan đến mật độ xương và rõ ràng hơn theo giới
Tuổi tăng 1 thì MDX giảm 0.002691, với chỉ số Adjusted R-squared: 0.2498
Việc 5. Xây dựng mô hình dự báo MĐX
5.1 Đọc dữ liệu ‘Oxteo-data.xlsx’ và gọi đối tượng là ‘fx’
library(readxl)## Warning: package 'readxl' was built under R version 4.4.1fx = as.data.frame(read_excelOsteo_data <- read_excel("C:/Dr Nhung/Khóa học phân tích dữ liệu trong NCKH bv 108 tháng 9.2024/Osteo-data.xlsx"))
dim(fx)## [1] 2216 17# Viết phương trình
summary(fx)## id sex age weight
## Min. : 1.0 Length:2216 Min. :57.00 Min. : 34.00
## 1st Qu.: 554.8 Class :character 1st Qu.:65.00 1st Qu.: 60.00
## Median :1108.5 Mode :character Median :70.00 Median : 69.00
## Mean :1108.5 Mean :70.89 Mean : 70.14
## 3rd Qu.:1662.2 3rd Qu.:76.00 3rd Qu.: 79.00
## Max. :2216.0 Max. :96.00 Max. :133.00
## NA's :53
## height prior_fx fnbmd smoking
## Min. :136.0 Min. :0.0000 Min. :0.2800 Min. :0.0000
## 1st Qu.:158.0 1st Qu.:0.0000 1st Qu.:0.7300 1st Qu.:0.0000
## Median :164.0 Median :0.0000 Median :0.8200 Median :0.0000
## Mean :164.9 Mean :0.1611 Mean :0.8287 Mean :0.4176
## 3rd Qu.:171.0 3rd Qu.:0.0000 3rd Qu.:0.9300 3rd Qu.:1.0000
## Max. :196.0 Max. :1.0000 Max. :1.5100 Max. :1.0000
## NA's :54 NA's :89 NA's :1
## parkinson rheum hypertension diabetes
## Min. :0.00000 Min. :0.00000 Min. :0.0000 Min. :0.000
## 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.000
## Median :0.00000 Median :0.00000 Median :1.0000 Median :0.000
## Mean :0.06498 Mean :0.03881 Mean :0.5063 Mean :0.111
## 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:1.0000 3rd Qu.:0.000
## Max. :1.00000 Max. :1.00000 Max. :1.0000 Max. :1.000
##
## copd cancer cvd falls_n
## Min. :0.000 Min. :0.00000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.000 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.0000
## Median :0.000 Median :0.00000 Median :0.0000 Median :0.0000
## Mean :0.111 Mean :0.08529 Mean :0.3872 Mean :0.2843
## 3rd Qu.:0.000 3rd Qu.:0.00000 3rd Qu.:1.0000 3rd Qu.:0.0000
## Max. :1.000 Max. :1.00000 Max. :1.0000 Max. :2.0000
##
## fx
## Min. :0.0000
## 1st Qu.:0.0000
## Median :0.0000
## Mean :0.2595
## 3rd Qu.:1.0000
## Max. :1.0000
## 5.2 Mô tả đặc điểm của dữ liệu theo giớ tính
5.3 Xây dựng mô hình dự báo MD0X tại cổ xương đùi
Phương pháp stepwise : Không khuyến cáo vì nhạy
fx.2 = na.omit(fx)
m.step = lm(fnbmd ~ age + sex + weight + height + fnbmd + prior_fx + falls_n + smoking + parkinson + rheum + hypertension + diabetes + copd + cancer + cvd, data = fx.2)## Warning in model.matrix.default(mt, mf, contrasts): the response appeared on
## the right-hand side and was dropped## Warning in model.matrix.default(mt, mf, contrasts): problem with term 5 in
## model.matrix: no columns are assignedstep = step(m.step)## Start: AIC=-9100.52
## fnbmd ~ age + sex + weight + height + fnbmd + prior_fx + falls_n +
## smoking + parkinson + rheum + hypertension + diabetes + copd +
## cancer + cvd## Warning in model.matrix.default(object, data = structure(list(fnbmd = c(1.08, :
## the response appeared on the right-hand side and was dropped## Warning in model.matrix.default(object, data = structure(list(fnbmd = c(1.08, :
## problem with term 5 in model.matrix: no columns are assigned##
## Step: AIC=-9100.52
## fnbmd ~ age + sex + weight + height + prior_fx + falls_n + smoking +
## parkinson + rheum + hypertension + diabetes + copd + cancer +
## cvd
##
## Df Sum of Sq RSS AIC
## - copd 1 0.0000 28.641 -9102.5
## - cvd 1 0.0029 28.643 -9102.3
## - hypertension 1 0.0037 28.644 -9102.2
## - cancer 1 0.0056 28.646 -9102.1
## - parkinson 1 0.0074 28.648 -9102.0
## - diabetes 1 0.0086 28.649 -9101.9
## - falls_n 1 0.0134 28.654 -9101.5
## <none> 28.641 -9100.5
## - rheum 1 0.0433 28.684 -9099.3
## - height 1 0.1607 28.801 -9090.7
## - prior_fx 1 0.3350 28.976 -9077.9
## - smoking 1 0.3807 29.021 -9074.5
## - sex 1 0.8234 29.464 -9042.4
## - age 1 2.1037 30.744 -8952.2
## - weight 1 4.9616 33.602 -8763.7
##
## Step: AIC=-9102.52
## fnbmd ~ age + sex + weight + height + prior_fx + falls_n + smoking +
## parkinson + rheum + hypertension + diabetes + cancer + cvd
##
## Df Sum of Sq RSS AIC
## - cvd 1 0.0029 28.643 -9104.3
## - hypertension 1 0.0037 28.644 -9104.2
## - cancer 1 0.0056 28.646 -9104.1
## - parkinson 1 0.0074 28.648 -9104.0
## - diabetes 1 0.0086 28.649 -9103.9
## - falls_n 1 0.0134 28.654 -9103.5
## <none> 28.641 -9102.5
## - rheum 1 0.0433 28.684 -9101.3
## - height 1 0.1607 28.801 -9092.7
## - prior_fx 1 0.3351 28.976 -9079.9
## - smoking 1 0.3810 29.022 -9076.5
## - sex 1 0.8236 29.464 -9044.4
## - age 1 2.1053 30.746 -8954.1
## - weight 1 4.9628 33.603 -8765.6
##
## Step: AIC=-9104.3
## fnbmd ~ age + sex + weight + height + prior_fx + falls_n + smoking +
## parkinson + rheum + hypertension + diabetes + cancer
##
## Df Sum of Sq RSS AIC
## - hypertension 1 0.0045 28.648 -9106.0
## - cancer 1 0.0057 28.649 -9105.9
## - parkinson 1 0.0074 28.651 -9105.8
## - diabetes 1 0.0076 28.651 -9105.7
## - falls_n 1 0.0134 28.657 -9105.3
## <none> 28.643 -9104.3
## - rheum 1 0.0426 28.686 -9103.2
## - height 1 0.1604 28.804 -9094.5
## - prior_fx 1 0.3346 28.978 -9081.7
## - smoking 1 0.3826 29.026 -9078.2
## - sex 1 0.8236 29.467 -9046.2
## - age 1 2.1048 30.748 -8955.9
## - weight 1 4.9620 33.605 -8767.4
##
## Step: AIC=-9105.97
## fnbmd ~ age + sex + weight + height + prior_fx + falls_n + smoking +
## parkinson + rheum + diabetes + cancer
##
## Df Sum of Sq RSS AIC
## - cancer 1 0.0057 28.654 -9107.5
## - diabetes 1 0.0066 28.655 -9107.5
## - parkinson 1 0.0077 28.656 -9107.4
## - falls_n 1 0.0130 28.661 -9107.0
## <none> 28.648 -9106.0
## - rheum 1 0.0428 28.691 -9104.8
## - height 1 0.1613 28.809 -9096.1
## - prior_fx 1 0.3334 28.981 -9083.4
## - smoking 1 0.3807 29.029 -9080.0
## - sex 1 0.8201 29.468 -9048.1
## - age 1 2.1035 30.751 -8957.7
## - weight 1 4.9596 33.608 -8769.3
##
## Step: AIC=-9107.55
## fnbmd ~ age + sex + weight + height + prior_fx + falls_n + smoking +
## parkinson + rheum + diabetes
##
## Df Sum of Sq RSS AIC
## - diabetes 1 0.0069 28.661 -9109.0
## - parkinson 1 0.0076 28.661 -9109.0
## - falls_n 1 0.0135 28.667 -9108.5
## <none> 28.654 -9107.5
## - rheum 1 0.0442 28.698 -9106.3
## - height 1 0.1628 28.817 -9097.5
## - prior_fx 1 0.3351 28.989 -9084.9
## - smoking 1 0.3854 29.039 -9081.2
## - sex 1 0.8219 29.476 -9049.6
## - age 1 2.0984 30.752 -8959.6
## - weight 1 4.9596 33.613 -8771.0
##
## Step: AIC=-9109.04
## fnbmd ~ age + sex + weight + height + prior_fx + falls_n + smoking +
## parkinson + rheum
##
## Df Sum of Sq RSS AIC
## - parkinson 1 0.0080 28.669 -9110.4
## - falls_n 1 0.0135 28.674 -9110.0
## <none> 28.661 -9109.0
## - rheum 1 0.0450 28.706 -9107.7
## - height 1 0.1614 28.822 -9099.1
## - prior_fx 1 0.3366 28.997 -9086.3
## - smoking 1 0.3844 29.045 -9082.8
## - sex 1 0.8262 29.487 -9050.8
## - age 1 2.1055 30.766 -8960.7
## - weight 1 4.9622 33.623 -8772.4
##
## Step: AIC=-9110.44
## fnbmd ~ age + sex + weight + height + prior_fx + falls_n + smoking +
## rheum
##
## Df Sum of Sq RSS AIC
## - falls_n 1 0.0132 28.682 -9111.5
## <none> 28.669 -9110.4
## - rheum 1 0.0454 28.714 -9109.1
## - height 1 0.1614 28.830 -9100.5
## - prior_fx 1 0.3345 29.003 -9087.8
## - smoking 1 0.3848 29.053 -9084.2
## - sex 1 0.8261 29.495 -9052.2
## - age 1 2.1123 30.781 -8961.7
## - weight 1 4.9672 33.636 -8773.5
##
## Step: AIC=-9111.47
## fnbmd ~ age + sex + weight + height + prior_fx + smoking + rheum
##
## Df Sum of Sq RSS AIC
## <none> 28.682 -9111.5
## - rheum 1 0.0467 28.729 -9110.0
## - height 1 0.1639 28.846 -9101.4
## - prior_fx 1 0.3383 29.020 -9088.6
## - smoking 1 0.3900 29.072 -9084.8
## - sex 1 0.8244 29.506 -9053.4
## - age 1 2.1100 30.792 -8962.9
## - weight 1 4.9608 33.643 -8775.15.3.2 Phương pháp BMA: Khuyến cáo dùng
library(BMA)## Warning: package 'BMA' was built under R version 4.4.1## Loading required package: survival## Loading required package: leaps## Warning: package 'leaps' was built under R version 4.4.1## Loading required package: robustbase## Warning: package 'robustbase' was built under R version 4.4.1##
## 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.4.1## Loading required package: rrcov## Warning: package 'rrcov' was built under R version 4.4.1## Scalable Robust Estimators with High Breakdown Point (version 1.7-6)xvars = fx.2[, c("age", "sex", "weight", "height", "prior_fx", "falls_n", "smoking", "parkinson", "rheum", "hypertension", "diabetes", "copd", "cancer", "cvd")]
m.bma = bicreg(xvars, fx.2$fnbmd, strict = FALSE, OR = 20)
summary(m.bma)##
## Call:
## bicreg(x = xvars, y = fx.2$fnbmd, strict = FALSE, OR = 20)
##
##
## 3 models were selected
## Best 3 models (cumulative posterior probability = 1 ):
##
## p!=0 EV SD model 1 model 2 model 3
## Intercept 100.0 0.626343 0.0977859 6.071e-01 6.098e-01 8.466e-01
## age 100.0 -0.004784 0.0003880 -4.765e-03 -4.778e-03 -4.995e-03
## sexMale 100.0 0.061554 0.0088666 6.021e-02 6.039e-02 7.690e-02
## weight 100.0 0.004329 0.0002379 4.306e-03 4.302e-03 4.603e-03
## height 92.1 0.001397 0.0005836 1.519e-03 1.503e-03 .
## prior_fx 100.0 -0.035389 0.0070901 -3.532e-02 -3.535e-02 -3.611e-02
## falls_n 0.0 0.000000 0.0000000 . . .
## smoking 100.0 -0.029091 0.0054107 -2.910e-02 -2.898e-02 -2.918e-02
## parkinson 0.0 0.000000 0.0000000 . . .
## rheum 10.0 0.002423 0.0083564 . 2.421e-02 .
## hypertension 0.0 0.000000 0.0000000 . . .
## diabetes 0.0 0.000000 0.0000000 . . .
## copd 0.0 0.000000 0.0000000 . . .
## cancer 0.0 0.000000 0.0000000 . . .
## cvd 0.0 0.000000 0.0000000 . . .
##
## nVar 6 7 5
## r2 0.434 0.435 0.431
## BIC -1.162e+03 -1.157e+03 -1.157e+03
## post prob 0.821 0.100 0.079post prob = sác xuất hậu định (%)
imageplot.bma(m.bma)
# 5.4 Kiểm tra giả định của mô hình
m.bmd = lm(fnbmd ~ age + sex + weight + height + prior_fx + smoking, data = fx.2)
summary(m.bmd)##
## Call:
## lm(formula = fnbmd ~ age + sex + weight + height + prior_fx +
## smoking, data = fx.2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.37425 -0.07568 -0.00663 0.07184 0.58750
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6070710 0.0758296 8.006 1.94e-15 ***
## age -0.0047646 0.0003834 -12.428 < 2e-16 ***
## sexMale 0.0602131 0.0077533 7.766 1.25e-14 ***
## weight 0.0043063 0.0002252 19.126 < 2e-16 ***
## height 0.0015189 0.0004328 3.509 0.000459 ***
## prior_fx -0.0353237 0.0070861 -4.985 6.70e-07 ***
## smoking -0.0290965 0.0054097 -5.379 8.34e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1166 on 2114 degrees of freedom
## Multiple R-squared: 0.4341, Adjusted R-squared: 0.4325
## F-statistic: 270.3 on 6 and 2114 DF, p-value: < 2.2e-16par(mfrow = c(2,2))
plot(m.bmd)
5.5 Viết công thức của mô hình tối ưu
summary(m.bmd)##
## Call:
## lm(formula = fnbmd ~ age + sex + weight + height + prior_fx +
## smoking, data = fx.2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.37425 -0.07568 -0.00663 0.07184 0.58750
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6070710 0.0758296 8.006 1.94e-15 ***
## age -0.0047646 0.0003834 -12.428 < 2e-16 ***
## sexMale 0.0602131 0.0077533 7.766 1.25e-14 ***
## weight 0.0043063 0.0002252 19.126 < 2e-16 ***
## height 0.0015189 0.0004328 3.509 0.000459 ***
## prior_fx -0.0353237 0.0070861 -4.985 6.70e-07 ***
## smoking -0.0290965 0.0054097 -5.379 8.34e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1166 on 2114 degrees of freedom
## Multiple R-squared: 0.4341, Adjusted R-squared: 0.4325
## F-statistic: 270.3 on 6 and 2114 DF, p-value: < 2.2e-16Phương trình:
MĐX = 0.607 - 0.0047646 x age + 0.0602131 x sex + … - 0.0290965 x smoking. Ex: BN nam, 70 tuổi, 70kg, 1,65 m, tiểu đường, ko hút thuốc…. áp vào phương trình (có thể làm trên R thường)