BMD and mortality Meta-analysis

(1) Analysis of individual projects

library(survival)
library(smoothHR)

## Loading required package: splines

library(table1)

## 
## Attaching package: 'table1'

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

(1.1) MrOS Study

Set-up dataset

mros = read.csv("C:\\Thach\\Research projects\\BMD and mortality\\Data\\MrOS_predictors.csv")
dim(mros)

## [1] 5994   35

head(mros)

##       ID B1FND entry_end age_0 AnyFx1st Hip Vert Proximal Distal entry_1stfx
## 1 BI0001 0.867 13.971253    67        0   0    0        0      0   13.971253
## 2 BI0002 0.740 14.636550    67        0   0    0        0      0   14.636550
## 3 BI0003 0.692 13.949350    72        1   0    0        0      1    8.421629
## 4 BI0004 0.878 18.825462    65        0   0    0        0      0   18.825462
## 5 BI0005 0.919  3.036277    78        0   0    0        0      0    3.036277
## 6 BI0006 0.886 13.149897    78        0   0    0        0      0   13.149897
##   age_1stfx AnyFx2nd entry_2ndfx age_2ndfx fx1st_2nd Death fx1st_end fx2nd_end
## 1      81.0        0   13.971253      81.0  0.000000     1  0.000000  0.000000
## 2      81.6        0   14.636550      81.6  0.000000     1  0.000000  0.000000
## 3      80.4        1   10.009582      82.0  1.587953     1  5.527721  3.939767
## 4      83.8        0   18.825462      83.8  0.000000     0  0.000000  0.000000
## 5      81.0        0    3.036277      81.0  0.000000     1  0.000000  0.000000
## 6      91.1        0   13.149897      91.1  0.000000     0  0.000000  0.000000
##   Age_death cvd diabetes copd cancer hypert parkinson ra activity smoke alcohol
## 1      81.0   0        0    0      0      0         1  0        1     0       0
## 2      81.6   1        0    0      0      1         1  1        1     1       0
## 3      85.9   1        0    0      0      0         0  0        1     1       1
## 4        NA   0        0    0      0      0         0  0        0     1       1
## 5      81.0   0        0    0      1      0         0  0        1     1       0
## 6        NA   0        0    1      0      0         0  1        0     1       0
##   alcohol_mod no_falls no_priorfx weight height      bmi
## 1           0        0          0   87.8 172.75 29.42107
## 2           0        0          0   99.4 176.60 31.87168
## 3           0        0          1   83.1 178.45 26.09564
## 4           0        0          0   97.4 184.65 28.56672
## 5           0        2          0   82.9 168.75 29.11166
## 6           0        1          0   76.3 163.75 28.45522

mros = subset(mros, select = c("ID", "B1FND", "entry_end", "Death", "age_0", "AnyFx1st", "cvd", "diabetes", "copd", "cancer", "hypert", "parkinson", "ra", "activity", "smoke", "alcohol", "bmi"))

(1.1.1) Baseline characteristics by survival status

table1(~ entry_end + B1FND + age_0 + as.factor(cvd) + as.factor(copd) + as.factor(diabetes) + as.factor(cancer) + as.factor(hypert) + as.factor(parkinson) + as.factor(ra) + as.factor(activity) + as.factor(smoke) + as.factor(alcohol) + bmi | Death, data = mros)

## Warning in table1.formula(~entry_end + B1FND + age_0 + as.factor(cvd) + : Terms
## to the right of '|' in formula 'x' define table columns and are expected to be
## factors with meaningful labels.

	0 (N=2382)	1 (N=3612)	Overall (N=5994)
entry_end
Mean (SD)	16.3 (3.78)	10.2 (4.73)	12.6 (5.30)
Median [Min, Max]	17.6 [0, 19.0]	10.5 [0.0383, 19.1]	13.9 [0, 19.1]
B1FND
Mean (SD)	0.796 (0.125)	0.776 (0.129)	0.784 (0.128)
Median [Min, Max]	0.784 [0.404, 1.49]	0.766 [0.273, 1.60]	0.774 [0.273, 1.60]
Missing	0 (0%)	1 (0.0%)	1 (0.0%)
age_0
Mean (SD)	70.8 (4.66)	75.5 (5.84)	73.7 (5.87)
Median [Min, Max]	70.0 [64.0, 91.0]	75.0 [65.0, 100]	73.0 [64.0, 100]
as.factor(cvd)
0	2071 (86.9%)	2672 (74.0%)	4743 (79.1%)
1	311 (13.1%)	940 (26.0%)	1251 (20.9%)
as.factor(copd)
0	2198 (92.3%)	3156 (87.4%)	5354 (89.3%)
1	184 (7.7%)	456 (12.6%)	640 (10.7%)
as.factor(diabetes)
0	2197 (92.2%)	3144 (87.0%)	5341 (89.1%)
1	185 (7.8%)	468 (13.0%)	653 (10.9%)
as.factor(cancer)
0	2046 (85.9%)	2856 (79.1%)	4902 (81.8%)
1	336 (14.1%)	756 (20.9%)	1092 (18.2%)
as.factor(hypert)
0	1507 (63.3%)	1906 (52.8%)	3413 (56.9%)
1	875 (36.7%)	1706 (47.2%)	2581 (43.1%)
as.factor(parkinson)
0	2082 (87.4%)	3123 (86.5%)	5205 (86.8%)
1	300 (12.6%)	489 (13.5%)	789 (13.2%)
as.factor(ra)
0	1353 (56.8%)	1794 (49.7%)	3147 (52.5%)
1	1029 (43.2%)	1818 (50.3%)	2847 (47.5%)
as.factor(activity)
0	646 (27.1%)	1355 (37.5%)	2001 (33.4%)
1	1736 (72.9%)	2257 (62.5%)	3993 (66.6%)
as.factor(smoke)
0	977 (41.0%)	1272 (35.2%)	2249 (37.5%)
1	1405 (59.0%)	2340 (64.8%)	3745 (62.5%)
as.factor(alcohol)
0	785 (33.0%)	1344 (37.2%)	2129 (35.5%)
1	1597 (67.0%)	2268 (62.8%)	3865 (64.5%)
bmi
Mean (SD)	27.4 (3.60)	27.4 (3.98)	27.4 (3.83)
Median [Min, Max]	27.0 [17.5, 50.7]	26.8 [17.2, 48.5]	26.9 [17.2, 50.7]
Missing	1 (0.0%)	1 (0.0%)	2 (0.0%)

(1.1.2) Age- adjusted model

mros$fnbmd = mros$B1FND/(-0.128)
mros.m1 <- coxph(Surv(entry_end, Death) ~ age_0 + fnbmd, data = mros, id = ID)
  summary(mros.m1)

## Call:
## coxph(formula = Surv(entry_end, Death) ~ age_0 + fnbmd, data = mros, 
##     id = ID)
## 
##   n= 5993, number of events= 3611 
##    (1 observation deleted due to missingness)
## 
##           coef exp(coef) se(coef)      z Pr(>|z|)    
## age_0 0.119908  1.127393 0.002994 40.049   <2e-16 ***
## fnbmd 0.020043  1.020245 0.017141  1.169    0.242    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## age_0     1.127     0.8870    1.1208     1.134
## fnbmd     1.020     0.9802    0.9865     1.055
## 
## Concordance= 0.682  (se = 0.005 )
## Likelihood ratio test= 1585  on 2 df,   p=<2e-16
## Wald test            = 1674  on 2 df,   p=<2e-16
## Score (logrank) test = 1762  on 2 df,   p=<2e-16

  zph.mros.m1 <- cox.zph(mros.m1)
  zph.mros.m1

##        chisq df     p
## age_0   3.91  1 0.048
## fnbmd   2.28  1 0.131
## GLOBAL  7.40  2 0.025

Examine the contribution of FNBMD to mortality risk

fit.1 = coxph(Surv(entry_end, Death) ~ age_0 + pspline(fnbmd), data = mros, x = TRUE)
hr.1 = smoothHR(data = mros, coxfit = fit.1)
plot(hr.1, predictor = "fnbmd", prob = 0.5, main = "MrOS Study - LogHR for FNBMD at median value")

(1.1.3) Multiple variable- adjusted model

mros.m2 <- coxph(Surv(entry_end, Death) ~ age_0 + fnbmd + cvd + diabetes + copd + cancer + hypert + ra + parkinson + activity + smoke + bmi, data = mros, id = ID)
  summary(mros.m2)

## Call:
## coxph(formula = Surv(entry_end, Death) ~ age_0 + fnbmd + cvd + 
##     diabetes + copd + cancer + hypert + ra + parkinson + activity + 
##     smoke + bmi, data = mros, id = ID)
## 
##   n= 5991, number of events= 3610 
##    (3 observations deleted due to missingness)
## 
##                coef exp(coef)  se(coef)      z Pr(>|z|)    
## age_0      0.116879  1.123984  0.003153 37.070  < 2e-16 ***
## fnbmd      0.059226  1.061014  0.018168  3.260 0.001115 ** 
## cvd        0.388356  1.474554  0.038903  9.983  < 2e-16 ***
## diabetes   0.373300  1.452520  0.051291  7.278 3.39e-13 ***
## copd       0.287483  1.333068  0.050593  5.682 1.33e-08 ***
## cancer     0.156760  1.169715  0.041442  3.783 0.000155 ***
## hypert     0.188418  1.207338  0.034474  5.465 4.62e-08 ***
## ra         0.063767  1.065844  0.034038  1.873 0.061015 .  
## parkinson  0.064600  1.066733  0.048895  1.321 0.186435    
## activity  -0.162974  0.849613  0.035137 -4.638 3.51e-06 ***
## smoke      0.238055  1.268778  0.035388  6.727 1.73e-11 ***
## bmi        0.012723  1.012805  0.005011  2.539 0.011111 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##           exp(coef) exp(-coef) lower .95 upper .95
## age_0        1.1240     0.8897    1.1171    1.1310
## fnbmd        1.0610     0.9425    1.0239    1.0995
## cvd          1.4746     0.6782    1.3663    1.5914
## diabetes     1.4525     0.6885    1.3136    1.6061
## copd         1.3331     0.7501    1.2072    1.4720
## cancer       1.1697     0.8549    1.0785    1.2687
## hypert       1.2073     0.8283    1.1285    1.2917
## ra           1.0658     0.9382    0.9971    1.1394
## parkinson    1.0667     0.9374    0.9693    1.1740
## activity     0.8496     1.1770    0.7931    0.9102
## smoke        1.2688     0.7882    1.1838    1.3599
## bmi          1.0128     0.9874    1.0029    1.0228
## 
## Concordance= 0.707  (se = 0.004 )
## Likelihood ratio test= 1985  on 12 df,   p=<2e-16
## Wald test            = 2050  on 12 df,   p=<2e-16
## Score (logrank) test = 2177  on 12 df,   p=<2e-16

  zph.mros.m2 <- cox.zph(mros.m2)
  zph.mros.m2

##             chisq df       p
## age_0      4.3705  1  0.0366
## fnbmd      2.4228  1  0.1196
## cvd        5.2895  1  0.0215
## diabetes   4.8684  1  0.0274
## copd       3.4565  1  0.0630
## cancer     2.7430  1  0.0977
## hypert     4.9692  1  0.0258
## ra         0.0014  1  0.9701
## parkinson  0.0809  1  0.7761
## activity   9.7407  1  0.0018
## smoke      3.3777  1  0.0661
## bmi        1.4601  1  0.2269
## GLOBAL    49.9770 12 1.4e-06

Examine the contribution of FNBMD to mortality risk

fit.2 = coxph(Surv(entry_end, Death) ~ age_0 + pspline(fnbmd)  + cvd + diabetes + copd + cancer + hypert + ra + parkinson + activity + smoke + bmi, data = mros, x = TRUE)
hr.2 = smoothHR(data = mros, coxfit = fit.2)

## Warning: a global p-value of  1.062086e-05 was obtained when testing the proportional hazards assumption

plot(hr.2, predictor = "fnbmd", prob = 0.5, main = "MrOS Study - Log HR for FNBMD at median value")

(1.2) SOF Study

Set-up dataset

sof = read.csv("C:\\Thach\\Research projects\\BMD and mortality\\Data\\SOF_death2.csv")
dim(sof)

## [1] 9408   76

head(sof)

##   ID Visit AGE WGHT  HGHT      BMI DAYS   TLD   FND      Tscore HPDAYS DEATH
## 1  1     2  69 67.3 150.5 29.71270  740 0.740 0.656 -1.67826087    740     0
## 2  3     2  75 72.7 151.0 31.88457  657 0.925 0.654 -1.69565217    657     0
## 3  4     2  67 60.6 155.1 25.19121  771 1.005 0.842 -0.06086957    771     0
## 4  5     2  84 44.7 157.0 18.13461  756 0.805 0.487 -3.14782609    756     0
## 5  6     2  71 84.1 161.6 32.20426  952 0.958 0.851  0.01739130    952     0
## 6  7     2  80 64.5 157.7 25.93560  706 0.585 0.560 -2.51304348    706     0
##   FOLALL MEG2 FALL NFALL ANYTOV    CA PHYS SMOKE DRWK30 CAL EST BISPH ECHF
## 1   6291   52    1     2      1    NA    0     1   1.25   1   0     0    0
## 2   6807   45    1     3      0 799.4    1     1   3.00   1   0     0    0
## 3   7213   45    1     1      0 889.2    1     0   0.25   1   0     0    0
## 4   5432   45    0     0      0 513.4    0     0   0.25   1   0     0    0
## 5   7165   46    0     0      0 235.3    1     0   0.00   1   1     0    0
## 6   1379   40    0     0      1    NA   NA     0   0.25  NA  NA    NA    0
##   ENGHRT EANGIN MURM HYTEN CVD EDIAB DIABCL EHTHY ECOPD REL EPARK EALZH ENEUR
## 1      0      0    0     0   0     0      0     0     0   0     0     0     0
## 2      0      0    0     0   0     0      0     0     1   0     0     0     0
## 3      0      0    0     0   0     0      0     0     0   0     0     0     0
## 4      0      0    0     0   0     0      0     0     0   1     0     0     0
## 5      0      0    0     1   0     0      0     0     0   0     1     0     0
## 6      0      0    0     0   0     0      0     0     0  NA     0    NA    NA
##   EDEPR ERA ECANCR ANYI ANYF HIPI HIPF PELI PELF VERTI VERTF HUMI HUMF RIBI
## 1     0   0      0    0 6291    0 6291    0 6291     0  4461    0 6291    0
## 2     0   0      1    0 6807    0 6807    0 6807     0  4453    0 6807    0
## 3     0   0      1   NA   NA    0 7213    0 7213     0  4595    0 7213    0
## 4     0   0      0    0 5432    0 5432    0 5432    NA  4610    0 5432    0
## 5     0   0      0    0 7165    0 7165    0 7165     0  4670    0 7165    0
## 6    NA  NA     NA    0 1379    0 1379    0 1379    NA  4916    0 1379    0
##   RIBF CLAVI CLAVF LLEGI LLEGF ELBI ELBF KNEEI KNEEF ULEGI ULEGF ANKI ANKF
## 1 6291     0  6291     0  6291    0 6291     0  6291     0  6291    0 6291
## 2 6807     0  6807     0  6807    0 6807     0  6807     0  6807    0 6807
## 3 7213     0  7213     0  7213    0 7213     0  7213     0  7213    0 7213
## 4 5432     0  5432     0  5432    0 5432     0  5432     0  5432    0 5432
## 5 7165     0  7165     0  7165    0 7165     0  7165     0  7165    0 7165
## 6 1379     0  1379     0  1379    0 1379     0  1379     0  1379    0 1379
##   HANDI HANDF FOOTI FOOTF WRSTI WRSTF VTTRI VTTRF Fx_TYPE Fx_FOLL Fx_v1v2
## 1     0  6291     0  6291     0  6291     0  6291       0    6291       0
## 2     0  6807     0  6807     0  6807     0  6807       0    6807       0
## 3     0  7213     0  7213     0  7213    NA    NA       0      NA       1
## 4     0  5432     0  5432     0  5432     0  5432       0    5432       0
## 5     0  7165     0  7165     0  7165     0  7165       0    7165       0
## 6     0  1379     0  1379     0  1379     0  1379       0    1379       0

sof = subset(sof, select = c("ID", "AGE", "BMI", "FND", "Tscore", "FOLALL", "DEATH", "ECHF", "EANGIN", "HYTEN", "CVD", "EDIAB", "ECOPD", "EPARK", "EALZH", "ENEUR", "EDEPR", "ERA", "ECANCR", "PHYS", "SMOKE", "DRWK30"))

sof$fu = sof$FOLALL/365.25
sd(sof$FND, na.rm = T)

## [1] 0.1105241

sof$fnbmd = sof$FND/(-0.11)

(1.2.1) Baseline characteristics by survival status

table1(~ fu + FND + fnbmd + Tscore + AGE + BMI + as.factor(CVD) + as.factor(ECHF) + as.factor(EANGIN) + as.factor(HYTEN) + as.factor(ECOPD) + as.factor(EPARK) + as.factor(ENEUR) + as.factor(EDEPR) + as.factor(ERA) + as.factor(ECANCR) + as.factor(EALZH) + as.factor(PHYS) + as.factor(SMOKE) + DRWK30 | DEATH, data = sof)

## Warning in table1.formula(~fu + FND + fnbmd + Tscore + AGE + BMI +
## as.factor(CVD) + : Terms to the right of '|' in formula 'x' define table
## columns and are expected to be factors with meaningful labels.

	0 (N=3743)	1 (N=5665)	Overall (N=9408)
fu
Mean (SD)	16.6 (4.54)	10.7 (5.16)	13.0 (5.72)
Median [Min, Max]	18.2 [0.107, 21.1]	11.0 [0.0329, 20.8]	13.8 [0.0329, 21.1]
Missing	70 (1.9%)	0 (0%)	70 (0.7%)
FND
Mean (SD)	0.662 (0.108)	0.640 (0.111)	0.649 (0.111)
Median [Min, Max]	0.652 [0.302, 1.26]	0.631 [0.277, 1.32]	0.639 [0.277, 1.32]
Missing	572 (15.3%)	876 (15.5%)	1448 (15.4%)
fnbmd
Mean (SD)	-6.02 (0.986)	-5.82 (1.01)	-5.90 (1.00)
Median [Min, Max]	-5.93 [-11.5, -2.75]	-5.74 [-12.0, -2.52]	-5.81 [-12.0, -2.52]
Missing	572 (15.3%)	876 (15.5%)	1448 (15.4%)
Tscore
Mean (SD)	-1.63 (0.943)	-1.82 (0.965)	-1.74 (0.961)
Median [Min, Max]	-1.71 [-4.76, 3.61]	-1.90 [-4.97, 4.13]	-1.83 [-4.97, 4.13]
Missing	572 (15.3%)	876 (15.5%)	1448 (15.4%)
AGE
Mean (SD)	71.6 (3.94)	74.9 (5.38)	73.6 (5.12)
Median [Min, Max]	71.0 [67.0, 89.0]	74.0 [67.0, 90.0]	72.0 [67.0, 90.0]
Missing	71 (1.9%)	47 (0.8%)	118 (1.3%)
BMI
Mean (SD)	26.2 (4.25)	26.2 (4.54)	26.2 (4.43)
Median [Min, Max]	25.6 [17.1, 42.0]	25.6 [16.9, 42.2]	25.6 [16.9, 42.2]
Missing	589 (15.7%)	823 (14.5%)	1412 (15.0%)
as.factor(CVD)
0	3140 (83.9%)	4199 (74.1%)	7339 (78.0%)
1	499 (13.3%)	1081 (19.1%)	1580 (16.8%)
2	91 (2.4%)	280 (4.9%)	371 (3.9%)
3	12 (0.3%)	84 (1.5%)	96 (1.0%)
4	1 (0.0%)	21 (0.4%)	22 (0.2%)
as.factor(ECHF)
0	3713 (99.2%)	5463 (96.4%)	9176 (97.5%)
1	30 (0.8%)	202 (3.6%)	232 (2.5%)
as.factor(EANGIN)
0	3483 (93.1%)	4962 (87.6%)	8445 (89.8%)
1	260 (6.9%)	703 (12.4%)	963 (10.2%)
as.factor(HYTEN)
0	2652 (70.9%)	3157 (55.7%)	5809 (61.7%)
1	1091 (29.1%)	2508 (44.3%)	3599 (38.3%)
as.factor(ECOPD)
0	3323 (88.8%)	4866 (85.9%)	8189 (87.0%)
1	225 (6.0%)	599 (10.6%)	824 (8.8%)
Missing	195 (5.2%)	200 (3.5%)	395 (4.2%)
as.factor(EPARK)
0	3729 (99.6%)	5622 (99.2%)	9351 (99.4%)
1	14 (0.4%)	43 (0.8%)	57 (0.6%)
as.factor(ENEUR)
0	3372 (90.1%)	4577 (80.8%)	7949 (84.5%)
1	38 (1.0%)	93 (1.6%)	131 (1.4%)
Missing	333 (8.9%)	995 (17.6%)	1328 (14.1%)
as.factor(EDEPR)
0	3161 (84.5%)	4207 (74.3%)	7368 (78.3%)
1	249 (6.7%)	465 (8.2%)	714 (7.6%)
Missing	333 (8.9%)	993 (17.5%)	1326 (14.1%)
as.factor(ERA)
0	3237 (86.5%)	4333 (76.5%)	7570 (80.5%)
1	172 (4.6%)	333 (5.9%)	505 (5.4%)
Missing	334 (8.9%)	999 (17.6%)	1333 (14.2%)
as.factor(ECANCR)
0	2759 (73.7%)	3588 (63.3%)	6347 (67.5%)
1	651 (17.4%)	1080 (19.1%)	1731 (18.4%)
Missing	333 (8.9%)	997 (17.6%)	1330 (14.1%)
as.factor(EALZH)
0	3399 (90.8%)	4595 (81.1%)	7994 (85.0%)
1	11 (0.3%)	78 (1.4%)	89 (0.9%)
Missing	333 (8.9%)	992 (17.5%)	1325 (14.1%)
as.factor(PHYS)
0	905 (24.2%)	1948 (34.4%)	2853 (30.3%)
1	2504 (66.9%)	2727 (48.1%)	5231 (55.6%)
Missing	334 (8.9%)	990 (17.5%)	1324 (14.1%)
as.factor(SMOKE)
0	2337 (62.4%)	3359 (59.3%)	5696 (60.5%)
1	1125 (30.1%)	1653 (29.2%)	2778 (29.5%)
2	281 (7.5%)	653 (11.5%)	934 (9.9%)
DRWK30
Mean (SD)	1.87 (3.70)	1.71 (3.75)	1.78 (3.73)
Median [Min, Max]	0.250 [0, 42.0]	0.250 [0, 42.0]	0.250 [0, 42.0]

(1.2.2) Age- adjusted model

sof.m1 <- coxph(Surv(fu, DEATH) ~ AGE + fnbmd, data = sof, id = ID)
  summary(sof.m1)

## Call:
## coxph(formula = Surv(fu, DEATH) ~ AGE + fnbmd, data = sof, id = ID)
## 
##   n= 7927, number of events= 4758 
##    (1481 observations deleted due to missingness)
## 
##           coef exp(coef) se(coef)      z Pr(>|z|)    
## AGE   0.118429  1.125727 0.002956 40.062  < 2e-16 ***
## fnbmd 0.056447  1.058071 0.015435  3.657 0.000255 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## AGE       1.126     0.8883     1.119     1.132
## fnbmd     1.058     0.9451     1.027     1.091
## 
## Concordance= 0.658  (se = 0.004 )
## Likelihood ratio test= 1611  on 2 df,   p=<2e-16
## Wald test            = 1804  on 2 df,   p=<2e-16
## Score (logrank) test = 1896  on 2 df,   p=<2e-16

  zph.sof.m1 <- cox.zph(sof.m1)
  zph.sof.m1

##        chisq df       p
## AGE    17.84  1 2.4e-05
## fnbmd   6.55  1   0.011
## GLOBAL 31.10  2 1.8e-07

Examine the contribution of FNBMD to mortality risk

fit.1 = coxph(Surv(fu, DEATH) ~ AGE + pspline(fnbmd), data = sof, x = TRUE)
hr.1 = smoothHR(data = sof, coxfit = fit.1)

## Warning: a global p-value of  3.378121e-06 was obtained when testing the proportional hazards assumption

plot(hr.1, predictor = "fnbmd", prob = 0.5, main = "SOF Study - LogHR for FNBMD at median value")

(1.2.3) Multiple variable- adjusted model

sof.m2 <- coxph(Surv(fu, DEATH) ~ AGE + fnbmd + BMI + CVD + EDIAB + ECOPD + EPARK + ENEUR + EDEPR + ERA + ECANCR + PHYS + SMOKE + DRWK30, data = sof, id = ID)
  summary(sof.m2)

## Call:
## coxph(formula = Surv(fu, DEATH) ~ AGE + fnbmd + BMI + CVD + EDIAB + 
##     ECOPD + EPARK + ENEUR + EDEPR + ERA + ECANCR + PHYS + SMOKE + 
##     DRWK30, data = sof, id = ID)
## 
##   n= 6813, number of events= 3940 
##    (2595 observations deleted due to missingness)
## 
##             coef exp(coef)  se(coef)      z Pr(>|z|)    
## AGE     0.119521  1.126957  0.003468 34.467  < 2e-16 ***
## fnbmd   0.063700  1.065773  0.018078  3.524 0.000426 ***
## BMI     0.007168  1.007194  0.004123  1.738 0.082125 .  
## CVD     0.212366  1.236600  0.025275  8.402  < 2e-16 ***
## EDIAB   0.592158  1.807886  0.060137  9.847  < 2e-16 ***
## ECOPD   0.262548  1.300238  0.053688  4.890 1.01e-06 ***
## EPARK   0.256515  1.292418  0.201207  1.275 0.202352    
## ENEUR   0.404095  1.497946  0.119968  3.368 0.000756 ***
## EDEPR   0.194627  1.214857  0.055668  3.496 0.000472 ***
## ERA     0.157103  1.170116  0.061765  2.544 0.010972 *  
## ECANCR  0.227291  1.255195  0.037898  5.997 2.01e-09 ***
## PHYS   -0.319798  0.726296  0.034454 -9.282  < 2e-16 ***
## SMOKE   0.263677  1.301708  0.025108 10.502  < 2e-16 ***
## DRWK30 -0.007415  0.992613  0.004662 -1.591 0.111706    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##        exp(coef) exp(-coef) lower .95 upper .95
## AGE       1.1270     0.8873    1.1193     1.135
## fnbmd     1.0658     0.9383    1.0287     1.104
## BMI       1.0072     0.9929    0.9991     1.015
## CVD       1.2366     0.8087    1.1768     1.299
## EDIAB     1.8079     0.5531    1.6069     2.034
## ECOPD     1.3002     0.7691    1.1704     1.445
## EPARK     1.2924     0.7737    0.8712     1.917
## ENEUR     1.4979     0.6676    1.1841     1.895
## EDEPR     1.2149     0.8231    1.0893     1.355
## ERA       1.1701     0.8546    1.0367     1.321
## ECANCR    1.2552     0.7967    1.1653     1.352
## PHYS      0.7263     1.3768    0.6789     0.777
## SMOKE     1.3017     0.7682    1.2392     1.367
## DRWK30    0.9926     1.0074    0.9836     1.002
## 
## Concordance= 0.694  (se = 0.004 )
## Likelihood ratio test= 1825  on 14 df,   p=<2e-16
## Wald test            = 2010  on 14 df,   p=<2e-16
## Score (logrank) test = 2110  on 14 df,   p=<2e-16

  zph.sof.m2 <- cox.zph(sof.m2)
  zph.sof.m2

##         chisq df       p
## AGE     9.836  1  0.0017
## fnbmd   3.003  1  0.0831
## BMI     0.485  1  0.4862
## CVD     1.289  1  0.2562
## EDIAB   3.945  1  0.0470
## ECOPD   3.068  1  0.0799
## EPARK   1.831  1  0.1761
## ENEUR   0.200  1  0.6544
## EDEPR   3.303  1  0.0692
## ERA     0.102  1  0.7499
## ECANCR  5.372  1  0.0205
## PHYS    4.427  1  0.0354
## SMOKE   2.332  1  0.1267
## DRWK30  0.749  1  0.3867
## GLOBAL 42.700 14 9.6e-05

Examine the contribution of FNBMD to mortality risk

fit.2 = coxph(Surv(fu, DEATH) ~ AGE + pspline(fnbmd) + BMI + CVD + EDIAB + ECOPD + EPARK + ENEUR + EDEPR + ERA + ECANCR + PHYS + SMOKE + DRWK30, data = sof, x = TRUE)
hr.2 = smoothHR(data = sof, coxfit = fit.2)

## Warning: a global p-value of  0.0003695908 was obtained when testing the proportional hazards assumption

plot(hr.2, predictor = "fnbmd", prob = 0.5, main = "SOF Study - Log HR for FNBMD at median value")

(1.3) DOES Study

(1.3.1) DOES Study - Men

Set-up dataset

does = read.csv("C:\\Thach\\Research projects\\BMD and mortality\\Data\\DOES_BMD_death.csv")
dim(does)

## [1] 3752   24

head(does)

##   BMI_0 age_0 neuro_0 rheum_0 hypert_0 diabetes_0 respi_0 cancer_0 cvd_0
## 1 27.34 73.07       0       0        0          1       0        0     1
## 2 26.13 67.39       0       0        1          0       0        0     0
## 3 27.34 68.47       0       0        1          0       0        0     1
## 4 24.06 62.10       1       0        1          0       0        0     0
## 5 24.06 60.92       0       0        0          0       0        0     0
## 6 26.13 75.81       0       0        1          0       0        0     1
##   prior_fx death Fx1st_any FNBMD_T StudyID sex date_entry Fx1st_date falls_n
## 1        0     1         0  0.3250       3   0  15JUN1989                  0
## 2        0     1         1 -0.2833       8   1  17APR1989  30JAN2009       0
## 3        0     1         0 -0.2250       9   0  12JUN1990                  0
## 4        0     1         0 -1.3417      10   1  04JUN1990                  0
## 5        0     1         0 -1.9083      23   0  08AUG1989                  0
## 6        0     1         0 -2.1417      24   1  03MAY1989                  0
##         doa last_visit   fnbmd   fu_death FNBMD_V1 FNBMD_raw
## 1 30JUN2018  01JAN1990 -0.3250  0.5475702    1.079 0.9745250
## 2 30JUN2018  31JUL2011  0.2833 22.2861054    0.966 0.8260040
## 3 30JUN2018  01JUL1995  0.2250  5.0513347    1.013 0.8991750
## 4 30JUN2018  30DEC2016  1.3417 26.5735797    0.839 0.6989960
## 5 30JUN2018  21JUL2010  1.9083 20.9500342    0.811 0.6685629
## 6 30JUN2018  16AUG2001  2.1417 12.2874743    0.743 0.6029960

does.m = subset(does, sex == 0, select = c("StudyID", "fu_death", "death", "age_0", "BMI_0", "neuro_0", "rheum_0", "diabetes_0", "respi_0", "cvd_0", "hypert_0", "cancer_0", "FNBMD_T", "fnbmd", "FNBMD_raw", "FNBMD_V1"))

sd(does.m$FNBMD_raw)

## [1] 0.1674037

does.m$fnbmd_sd = does.m$FNBMD_raw/(-0.167)

(1.3.1.1) Baseline characteristics by survival status

table1(~ fu_death + FNBMD_raw + fnbmd + FNBMD_T + age_0 + BMI_0 + as.factor(cvd_0) + as.factor(diabetes_0) + as.factor(respi_0) + as.factor(neuro_0) + as.factor(cancer_0) + as.factor(hypert_0) + as.factor(rheum_0) | death, data = does.m)

## Warning in table1.formula(~fu_death + FNBMD_raw + fnbmd + FNBMD_T + age_0 + :
## Terms to the right of '|' in formula 'x' define table columns and are expected
## to be factors with meaningful labels.

	0 (N=616)	1 (N=745)	Overall (N=1361)
fu_death
Mean (SD)	17.3 (6.49)	11.0 (7.05)	13.8 (7.50)
Median [Min, Max]	15.8 [8.00, 29.2]	9.50 [0.0876, 28.2]	13.5 [0.0876, 29.2]
FNBMD_raw
Mean (SD)	0.816 (0.152)	0.781 (0.178)	0.797 (0.167)
Median [Min, Max]	0.818 [0.349, 1.36]	0.771 [0.103, 1.43]	0.795 [0.103, 1.43]
fnbmd
Mean (SD)	0.831 (1.11)	1.08 (1.30)	0.970 (1.22)
Median [Min, Max]	0.817 [-3.12, 4.24]	1.16 [-3.62, 6.03]	0.983 [-3.62, 6.03]
FNBMD_T
Mean (SD)	-0.831 (1.11)	-1.08 (1.30)	-0.970 (1.22)
Median [Min, Max]	-0.817 [-4.24, 3.12]	-1.16 [-6.03, 3.62]	-0.983 [-6.03, 3.62]
age_0
Mean (SD)	67.4 (4.76)	71.4 (6.41)	69.6 (6.05)
Median [Min, Max]	66.5 [59.6, 87.1]	70.0 [59.8, 91.9]	68.2 [59.6, 91.9]
BMI_0
Mean (SD)	27.7 (3.83)	27.4 (2.30)	27.5 (3.09)
Median [Min, Max]	27.6 [17.5, 43.4]	27.3 [15.6, 41.0]	27.3 [15.6, 43.4]
as.factor(cvd_0)
0	424 (68.8%)	403 (54.1%)	827 (60.8%)
1	192 (31.2%)	342 (45.9%)	534 (39.2%)
as.factor(diabetes_0)
0	522 (84.7%)	651 (87.4%)	1173 (86.2%)
1	94 (15.3%)	94 (12.6%)	188 (13.8%)
as.factor(respi_0)
0	545 (88.5%)	646 (86.7%)	1191 (87.5%)
1	71 (11.5%)	99 (13.3%)	170 (12.5%)
as.factor(neuro_0)
0	586 (95.1%)	695 (93.3%)	1281 (94.1%)
1	30 (4.9%)	50 (6.7%)	80 (5.9%)
as.factor(cancer_0)
0	561 (91.1%)	677 (90.9%)	1238 (91.0%)
1	55 (8.9%)	68 (9.1%)	123 (9.0%)
as.factor(hypert_0)
0	335 (54.4%)	416 (55.8%)	751 (55.2%)
1	281 (45.6%)	329 (44.2%)	610 (44.8%)
as.factor(rheum_0)
0	610 (99.0%)	724 (97.2%)	1334 (98.0%)
1	6 (1.0%)	21 (2.8%)	27 (2.0%)

(1.3.1.2) Age- adjusted model

does.m1 <- coxph(Surv(fu_death, death) ~ age_0 + fnbmd_sd, data = does.m, id = StudyID)
  summary(does.m1)

## Call:
## coxph(formula = Surv(fu_death, death) ~ age_0 + fnbmd_sd, data = does.m, 
##     id = StudyID)
## 
##   n= 1361, number of events= 745 
## 
##              coef exp(coef) se(coef) robust se     z Pr(>|z|)    
## age_0    0.097982  1.102943 0.005949  0.007500 13.06   <2e-16 ***
## fnbmd_sd 0.051434  1.052780 0.040404  0.051954  0.99    0.322    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## age_0        1.103     0.9067    1.0868     1.119
## fnbmd_sd     1.053     0.9499    0.9509     1.166
## 
## Concordance= 0.689  (se = 0.012 )
## Likelihood ratio test= 270.4  on 2 df,   p=<2e-16
## Wald test            = 199.5  on 2 df,   p=<2e-16
## Score (logrank) test = 319.5  on 2 df,   p=<2e-16,   Robust = 173.1  p=<2e-16
## 
##   (Note: the likelihood ratio and score tests assume independence of
##      observations within a cluster, the Wald and robust score tests do not).

  zph.does.m1 <- cox.zph(does.m1)
  zph.does.m1

##          chisq df       p
## age_0     11.2  1 0.00081
## fnbmd_sd  19.6  1 9.4e-06
## GLOBAL    26.1  2 2.2e-06

Examine the contribution of FNBMD to mortality risk

fit.1 = coxph(Surv(fu_death, death) ~ age_0 + pspline(fnbmd_sd), data = does.m, x = TRUE)
hr.1 = smoothHR(data = does.m, coxfit = fit.1)

## Warning: a global p-value of  0.0003276478 was obtained when testing the proportional hazards assumption

plot(hr.1, predictor = "fnbmd_sd", prob = 0.5, main = "DOES Men - LogHR for FNBMD at median value")

(1.3.1.3) Multiple variable- adjusted model

does.m2 <- coxph(Surv(fu_death, death) ~ age_0 + fnbmd_sd + BMI_0 + cvd_0 + diabetes_0 + respi_0 + cancer_0 + rheum_0 + hypert_0 , data = does.m, id = StudyID)
  summary(does.m2)

## Call:
## coxph(formula = Surv(fu_death, death) ~ age_0 + fnbmd_sd + BMI_0 + 
##     cvd_0 + diabetes_0 + respi_0 + cancer_0 + rheum_0 + hypert_0, 
##     data = does.m, id = StudyID)
## 
##   n= 1361, number of events= 745 
## 
##                 coef exp(coef)  se(coef) robust se      z Pr(>|z|)    
## age_0       0.098404  1.103409  0.006068  0.007752 12.694  < 2e-16 ***
## fnbmd_sd    0.099376  1.104482  0.042408  0.054420  1.826 0.067837 .  
## BMI_0       0.047735  1.048893  0.011303  0.013344  3.577 0.000347 ***
## cvd_0       0.203534  1.225727  0.075155  0.083604  2.434 0.014913 *  
## diabetes_0 -0.053323  0.948074  0.113194  0.129969 -0.410 0.681607    
## respi_0     0.091454  1.095766  0.108721  0.149186  0.613 0.539864    
## cancer_0    0.087792  1.091761  0.128362  0.133162  0.659 0.509709    
## rheum_0     0.277581  1.319933  0.222801  0.213289  1.301 0.193112    
## hypert_0   -0.051033  0.950247  0.077421  0.087848 -0.581 0.561286    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## age_0         1.1034     0.9063    1.0868     1.120
## fnbmd_sd      1.1045     0.9054    0.9927     1.229
## BMI_0         1.0489     0.9534    1.0218     1.077
## cvd_0         1.2257     0.8158    1.0405     1.444
## diabetes_0    0.9481     1.0548    0.7349     1.223
## respi_0       1.0958     0.9126    0.8180     1.468
## cancer_0      1.0918     0.9160    0.8410     1.417
## rheum_0       1.3199     0.7576    0.8690     2.005
## hypert_0      0.9502     1.0524    0.7999     1.129
## 
## Concordance= 0.692  (se = 0.012 )
## Likelihood ratio test= 298.5  on 9 df,   p=<2e-16
## Wald test            = 224.1  on 9 df,   p=<2e-16
## Score (logrank) test = 347.8  on 9 df,   p=<2e-16,   Robust = 189.9  p=<2e-16
## 
##   (Note: the likelihood ratio and score tests assume independence of
##      observations within a cluster, the Wald and robust score tests do not).

  zph.does.m2 <- cox.zph(does.m2)
  zph.does.m2

##             chisq df       p
## age_0      10.574  1  0.0011
## fnbmd_sd   15.789  1 7.1e-05
## BMI_0       0.384  1  0.5357
## cvd_0       4.012  1  0.0452
## diabetes_0  3.281  1  0.0701
## respi_0     4.113  1  0.0426
## cancer_0    7.719  1  0.0055
## rheum_0     4.186  1  0.0408
## hypert_0   22.275  1 2.4e-06
## GLOBAL     54.507  9 1.5e-08

Examine the contribution of FNBMD to mortality risk

fit.2 = coxph(Surv(fu_death, death) ~ age_0 + fnbmd_sd + BMI_0 + cvd_0 + diabetes_0 + respi_0 + cancer_0 + rheum_0 + hypert_0 , data = does.m, x = TRUE)

hr.2 = smoothHR(data = does.m, coxfit = fit.2)

## Warning: a global p-value of  1.512141e-08 was obtained when testing the proportional hazards assumption

plot(hr.2, predictor = "fnbmd", prob = 0.5, main = "DOES Men - Log HR for FNBMD at median value")

(1.3.2) DOES Study - Women

Set-up dataset

does.w = subset(does, sex == 1, select = c("StudyID", "fu_death", "death", "age_0", "BMI_0", "neuro_0", "rheum_0", "diabetes_0", "respi_0", "cvd_0", "hypert_0", "cancer_0", "FNBMD_T", "fnbmd", "FNBMD_raw", "FNBMD_V1"))

sd(does.w$FNBMD_raw)

## [1] 0.1379924

does.w$fnbmd_sd = does.w$FNBMD_raw/(-0.138)

(1.3.2.1) Baseline characteristics by survival status

table1(~ fu_death + FNBMD_raw + fnbmd + FNBMD_T + age_0 + BMI_0 + as.factor(cvd_0) + as.factor(diabetes_0) + as.factor(respi_0) + as.factor(neuro_0) + as.factor(cancer_0) + as.factor(hypert_0) + as.factor(rheum_0) | death, data = does.w)

## Warning in table1.formula(~fu_death + FNBMD_raw + fnbmd + FNBMD_T + age_0 + :
## Terms to the right of '|' in formula 'x' define table columns and are expected
## to be factors with meaningful labels.

	0 (N=1356)	1 (N=1035)	Overall (N=2391)
fu_death
Mean (SD)	17.4 (6.23)	12.8 (7.07)	15.4 (6.99)
Median [Min, Max]	16.0 [8.00, 29.5]	12.3 [0.0301, 28.9]	14.7 [0.0301, 29.5]
FNBMD_raw
Mean (SD)	0.699 (0.129)	0.629 (0.139)	0.669 (0.138)
Median [Min, Max]	0.688 [0.389, 1.15]	0.625 [0.139, 1.27]	0.663 [0.139, 1.27]
fnbmd
Mean (SD)	1.35 (1.08)	1.92 (1.16)	1.59 (1.15)
Median [Min, Max]	1.43 [-2.43, 3.93]	1.96 [-3.45, 6.01]	1.64 [-3.45, 6.01]
FNBMD_T
Mean (SD)	-1.35 (1.08)	-1.92 (1.16)	-1.59 (1.15)
Median [Min, Max]	-1.43 [-3.93, 2.43]	-1.96 [-6.01, 3.45]	-1.64 [-6.01, 3.45]
age_0
Mean (SD)	67.1 (5.22)	72.4 (7.52)	69.4 (6.83)
Median [Min, Max]	66.0 [59.7, 90.2]	71.5 [59.8, 94.1]	67.7 [59.7, 94.1]
BMI_0
Mean (SD)	27.1 (5.15)	26.3 (2.80)	26.8 (4.31)
Median [Min, Max]	26.6 [16.4, 53.7]	26.1 [15.4, 48.7]	26.1 [15.4, 53.7]
as.factor(cvd_0)
0	1004 (74.0%)	653 (63.1%)	1657 (69.3%)
1	352 (26.0%)	382 (36.9%)	734 (30.7%)
as.factor(diabetes_0)
0	1223 (90.2%)	903 (87.2%)	2126 (88.9%)
1	133 (9.8%)	132 (12.8%)	265 (11.1%)
as.factor(respi_0)
0	1198 (88.3%)	925 (89.4%)	2123 (88.8%)
1	158 (11.7%)	110 (10.6%)	268 (11.2%)
as.factor(neuro_0)
0	1256 (92.6%)	960 (92.8%)	2216 (92.7%)
1	100 (7.4%)	75 (7.2%)	175 (7.3%)
as.factor(cancer_0)
0	1225 (90.3%)	947 (91.5%)	2172 (90.8%)
1	131 (9.7%)	88 (8.5%)	219 (9.2%)
as.factor(hypert_0)
0	657 (48.5%)	483 (46.7%)	1140 (47.7%)
1	699 (51.5%)	552 (53.3%)	1251 (52.3%)
as.factor(rheum_0)
0	1312 (96.8%)	984 (95.1%)	2296 (96.0%)
1	44 (3.2%)	51 (4.9%)	95 (4.0%)

(1.3.2.2) Age- adjusted model

does.w1 <- coxph(Surv(fu_death, death) ~ age_0 + fnbmd_sd, data = does.w, id = StudyID)
  summary(does.w1)

## Call:
## coxph(formula = Surv(fu_death, death) ~ age_0 + fnbmd_sd, data = does.w, 
##     id = StudyID)
## 
##   n= 2391, number of events= 1035 
## 
##              coef exp(coef) se(coef) robust se      z Pr(>|z|)    
## age_0    0.107842  1.113872 0.004822  0.005823 18.522  < 2e-16 ***
## fnbmd_sd 0.133908  1.143288 0.038026  0.047723  2.806  0.00502 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## age_0        1.114     0.8978     1.101     1.127
## fnbmd_sd     1.143     0.8747     1.041     1.255
## 
## Concordance= 0.733  (se = 0.009 )
## Likelihood ratio test= 617.1  on 2 df,   p=<2e-16
## Wald test            = 435.7  on 2 df,   p=<2e-16
## Score (logrank) test = 751.8  on 2 df,   p=<2e-16,   Robust = 319.8  p=<2e-16
## 
##   (Note: the likelihood ratio and score tests assume independence of
##      observations within a cluster, the Wald and robust score tests do not).

  zph.does.w1 <- cox.zph(does.w1)
  zph.does.w1

##          chisq df       p
## age_0     41.3  1 1.3e-10
## fnbmd_sd  38.4  1 5.7e-10
## GLOBAL    56.0  2 7.0e-13

Examine the contribution of FNBMD to mortality risk

fit.1 = coxph(Surv(fu_death, death) ~ age_0 + pspline(fnbmd_sd), data = does.w, x = TRUE)
hr.1 = smoothHR(data = does.w, coxfit = fit.1)

## Warning: a global p-value of  7.19869e-09 was obtained when testing the proportional hazards assumption

plot(hr.1, predictor = "fnbmd_sd", prob = 0.5, main = "DOES Women - LogHR for FNBMD at median value")

(1.3.2.3) Multiple variable- adjusted model

does.w2 <- coxph(Surv(fu_death, death) ~ age_0 + fnbmd_sd + BMI_0 + cvd_0 + diabetes_0 + respi_0 + cancer_0 + rheum_0 + hypert_0 , data = does.w, id = StudyID)
  summary(does.w2)

## Call:
## coxph(formula = Surv(fu_death, death) ~ age_0 + fnbmd_sd + BMI_0 + 
##     cvd_0 + diabetes_0 + respi_0 + cancer_0 + rheum_0 + hypert_0, 
##     data = does.w, id = StudyID)
## 
##   n= 2391, number of events= 1035 
## 
##                 coef exp(coef)  se(coef) robust se      z Pr(>|z|)    
## age_0       0.109543  1.115768  0.004949  0.005880 18.630  < 2e-16 ***
## fnbmd_sd    0.173357  1.189290  0.040906  0.052540  3.300 0.000969 ***
## BMI_0       0.020226  1.020432  0.008402  0.008740  2.314 0.020650 *  
## cvd_0       0.034568  1.035173  0.065730  0.076738  0.450 0.652373    
## diabetes_0  0.440375  1.553289  0.096793  0.131423  3.351 0.000806 ***
## respi_0     0.020655  1.020870  0.101067  0.125403  0.165 0.869174    
## cancer_0    0.024937  1.025250  0.112391  0.122101  0.204 0.838174    
## rheum_0     0.373262  1.452465  0.144927  0.177641  2.101 0.035622 *  
## hypert_0   -0.037742  0.962962  0.063960  0.076024 -0.496 0.619582    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## age_0          1.116     0.8962    1.1030     1.129
## fnbmd_sd       1.189     0.8408    1.0729     1.318
## BMI_0          1.020     0.9800    1.0031     1.038
## cvd_0          1.035     0.9660    0.8906     1.203
## diabetes_0     1.553     0.6438    1.2006     2.010
## respi_0        1.021     0.9796    0.7984     1.305
## cancer_0       1.025     0.9754    0.8070     1.302
## rheum_0        1.452     0.6885    1.0254     2.057
## hypert_0       0.963     1.0385    0.8297     1.118
## 
## Concordance= 0.741  (se = 0.009 )
## Likelihood ratio test= 653.5  on 9 df,   p=<2e-16
## Wald test            = 466.8  on 9 df,   p=<2e-16
## Score (logrank) test = 790.4  on 9 df,   p=<2e-16,   Robust = 332.1  p=<2e-16
## 
##   (Note: the likelihood ratio and score tests assume independence of
##      observations within a cluster, the Wald and robust score tests do not).

  zph.does.w2 <- cox.zph(does.w2)
  zph.does.w2

##             chisq df       p
## age_0      40.299  1 2.2e-10
## fnbmd_sd   37.654  1 8.4e-10
## BMI_0       4.973  1 0.02575
## cvd_0       0.444  1 0.50519
## diabetes_0  3.711  1 0.05405
## respi_0     1.642  1 0.20000
## cancer_0    3.734  1 0.05333
## rheum_0    12.256  1 0.00046
## hypert_0   14.707  1 0.00013
## GLOBAL     90.192  9 1.5e-15

Examine the contribution of FNBMD to mortality risk

fit.2 = coxph(Surv(fu_death, death) ~ age_0 + fnbmd_sd + BMI_0 + cvd_0 + diabetes_0 + respi_0 + cancer_0 + rheum_0 + hypert_0 , data = does.w, x = TRUE)

hr.2 = smoothHR(data = does.w, coxfit = fit.2)

## Warning: a global p-value of  1.49014e-15 was obtained when testing the proportional hazards assumption

plot(hr.2, predictor = "fnbmd", prob = 0.5, main = "DOES Women - Log HR for FNBMD at median value")

BMD and mortality

Thach Tran

2023-07-04

BMD and mortality Meta-analysis

(1) Analysis of individual projects

(1.1) MrOS Study

(1.1.1) Baseline characteristics by survival status

(1.1.2) Age- adjusted model

(1.1.3) Multiple variable- adjusted model

(1.2) SOF Study

(1.2.1) Baseline characteristics by survival status

(1.2.2) Age- adjusted model

(1.2.3) Multiple variable- adjusted model

(1.3) DOES Study

(1.3.1) DOES Study - Men

(1.3.1.1) Baseline characteristics by survival status

(1.3.1.2) Age- adjusted model

(1.3.1.3) Multiple variable- adjusted model

(1.3.2) DOES Study - Women

(1.3.2.1) Baseline characteristics by survival status

(1.3.2.2) Age- adjusted model

(1.3.2.3) Multiple variable- adjusted model