day5
Cúc
2024-09-11
library(readxl)
fx = as.data.frame(read_excel("//Users//nguyenthicuc//Downloads//Osteo-data.xlsx"))
dim(fx)## [1] 2216 17summary(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
## library(table1)##
## Attaching package: 'table1'## The following objects are masked from 'package:base':
##
## units, units<-table1(~ age + weight + height + fnbmd + as.factor(prior_fx) + as.factor(falls_n) + as.factor(smoking) + as.factor(parkinson) + as.factor(rheum) + as.factor(hypertension) + as.factor(diabetes) + as.factor(copd) + as.factor(cancer) + as.factor(cvd) + as.factor(fx) | sex, data = fx)| Female (N=1358) |
Male (N=858) |
Overall (N=2216) |
|
|---|---|---|---|
| age | |||
| Mean (SD) | 71.2 (7.59) | 70.4 (6.44) | 70.9 (7.17) |
| Median [Min, Max] | 70.0 [57.0, 96.0] | 69.0 [59.0, 92.0] | 70.0 [57.0, 96.0] |
| weight | |||
| Mean (SD) | 64.9 (12.5) | 78.2 (12.7) | 70.1 (14.2) |
| Median [Min, Max] | 64.0 [34.0, 115] | 78.0 [45.0, 133] | 69.0 [34.0, 133] |
| Missing | 42 (3.1%) | 11 (1.3%) | 53 (2.4%) |
| height | |||
| Mean (SD) | 160 (6.37) | 173 (6.86) | 165 (9.35) |
| Median [Min, Max] | 160 [136, 181] | 173 [151, 196] | 164 [136, 196] |
| Missing | 41 (3.0%) | 13 (1.5%) | 54 (2.4%) |
| fnbmd | |||
| Mean (SD) | 0.777 (0.132) | 0.909 (0.153) | 0.829 (0.155) |
| Median [Min, Max] | 0.770 [0.280, 1.31] | 0.900 [0.340, 1.51] | 0.820 [0.280, 1.51] |
| Missing | 57 (4.2%) | 32 (3.7%) | 89 (4.0%) |
| as.factor(prior_fx) | |||
| 0 | 1125 (82.8%) | 734 (85.5%) | 1859 (83.9%) |
| 1 | 233 (17.2%) | 124 (14.5%) | 357 (16.1%) |
| as.factor(falls_n) | |||
| 0 | 1063 (78.3%) | 671 (78.2%) | 1734 (78.2%) |
| 1 | 206 (15.2%) | 128 (14.9%) | 334 (15.1%) |
| 2 | 89 (6.6%) | 59 (6.9%) | 148 (6.7%) |
| as.factor(smoking) | |||
| 0 | 962 (70.8%) | 328 (38.2%) | 1290 (58.2%) |
| 1 | 395 (29.1%) | 530 (61.8%) | 925 (41.7%) |
| Missing | 1 (0.1%) | 0 (0%) | 1 (0.0%) |
| as.factor(parkinson) | |||
| 0 | 1268 (93.4%) | 804 (93.7%) | 2072 (93.5%) |
| 1 | 90 (6.6%) | 54 (6.3%) | 144 (6.5%) |
| as.factor(rheum) | |||
| 0 | 1306 (96.2%) | 824 (96.0%) | 2130 (96.1%) |
| 1 | 52 (3.8%) | 34 (4.0%) | 86 (3.9%) |
| as.factor(hypertension) | |||
| 0 | 695 (51.2%) | 399 (46.5%) | 1094 (49.4%) |
| 1 | 663 (48.8%) | 459 (53.5%) | 1122 (50.6%) |
| as.factor(diabetes) | |||
| 0 | 1213 (89.3%) | 757 (88.2%) | 1970 (88.9%) |
| 1 | 145 (10.7%) | 101 (11.8%) | 246 (11.1%) |
| as.factor(copd) | |||
| 0 | 1211 (89.2%) | 759 (88.5%) | 1970 (88.9%) |
| 1 | 147 (10.8%) | 99 (11.5%) | 246 (11.1%) |
| as.factor(cancer) | |||
| 0 | 1235 (90.9%) | 792 (92.3%) | 2027 (91.5%) |
| 1 | 123 (9.1%) | 66 (7.7%) | 189 (8.5%) |
| as.factor(cvd) | |||
| 0 | 843 (62.1%) | 515 (60.0%) | 1358 (61.3%) |
| 1 | 515 (37.9%) | 343 (40.0%) | 858 (38.7%) |
| as.factor(fx) | |||
| 0 | 932 (68.6%) | 709 (82.6%) | 1641 (74.1%) |
| 1 | 426 (31.4%) | 149 (17.4%) | 575 (25.9%) |
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.1summary(step)##
## Call:
## lm(formula = fnbmd ~ age + sex + weight + height + prior_fx +
## smoking + rheum, data = fx.2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.37479 -0.07541 -0.00655 0.06870 0.56423
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6097744 0.0757999 8.045 1.43e-15 ***
## age -0.0047778 0.0003832 -12.468 < 2e-16 ***
## sexMale 0.0603920 0.0077495 7.793 1.02e-14 ***
## weight 0.0043022 0.0002250 19.117 < 2e-16 ***
## height 0.0015035 0.0004327 3.475 0.000521 ***
## prior_fx -0.0353546 0.0070820 -4.992 6.46e-07 ***
## smoking -0.0289828 0.0054069 -5.360 9.21e-08 ***
## rheum 0.0242102 0.0130528 1.855 0.063765 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1165 on 2113 degrees of freedom
## Multiple R-squared: 0.435, Adjusted R-squared: 0.4332
## F-statistic: 232.4 on 7 and 2113 DF, p-value: < 2.2e-16library(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-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.079imageplot.bma(m.bma)
xvars = fx.2[, c("age", "sex", "weight", "height", "prior_fx", "falls_n", "smoking", "parkinson", "rheum", "hypertension", "diabetes", "copd", "cancer", "cvd")]
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)
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-16