library(table1)
## 
## Attaching package: 'table1'
## The following objects are masked from 'package:base':
## 
##     units, units<-
library(compareGroups)
library(survival)

1. MEN (MrOS data)

1.1 Data setup

men = read.csv("C:\\UTS\\Research projects\\Skeletal Age_BMD and Mortality\\MrOS data.csv")
head(men)
##       ID B1FND GIAGE1 FUCDTIME COMORBIDITIES prior_nt_fx FALLS prior_fx_
## 1 SD8108 0.357     72     1984             1           1     0         1
## 2 PI4554 0.404     70     1829             2           0     0         0
## 3 SD8293 0.463     69     6657             0           0     0         0
## 4 PO7118 0.469     77     4239             1           0     0         0
## 5 PO6818 0.475     77     1376             3           0     1         0
## 6 PI5172 0.478     78     4824             0           0     0         0
##   DADEAD2 FX_TYPE entry_end age_0 AnyFx1st Hip Vert Proximal Distal entry_1stfx
## 1       1       A  5.431896    72        0   0    0        0      0    5.431896
## 2       0       A  5.007529    70        0   0    0        0      0    5.007529
## 3       0       A 18.225873    69        0   0    0        0      0   18.225873
## 4       1       A 11.605749    77        0   0    0        0      0   11.605749
## 5       1       A  3.767283    77        0   0    0        0      0    3.767283
## 6       1       A 13.207392    78        0   0    0        0      0   13.207392
##   age_1stfx AnyFx2nd entry_2ndfx age_2ndfx fx1st_2nd Death fx1st_end fx2nd_end
## 1      77.4        0    5.431896      77.4         0     1         0         0
## 2      75.0        0    5.007529      75.0         0     0         0         0
## 3      87.2        0   18.225873      87.2         0     0         0         0
## 4      88.6        0   11.605749      88.6         0     1         0         0
## 5      80.8        0    3.767283      80.8         0     1         0         0
## 6      91.2        0   13.207392      91.2         0     1         0         0
##   Age_death cvd_angina cvd diabetes copd cancer hypert parkinson ra activity
## 1      77.4          1   0        0    0      1      0         0  0        0
## 2        NA          1   1        0    0      0      1         0  1        0
## 3        NA          0   0        0    0      0      0         0  0        1
## 4      88.6          0   0        0    1      0      0         0  0        0
## 5      80.8          1   1        0    0      0      1         0  1        1
## 6      91.2          0   0        0    0      0      0         0  1        1
##   smoke alcohol alcohol_mod group      FNT no_falls Fall no_priorfx priorfx
## 1     1       1           1     0 4.178844        0    0          1       1
## 2     1       1           0     0 3.783882        0    0          0       0
## 3     0       1           0     0 3.295253        0    0          0       0
## 4     1       1           0     0 3.242558        0    0          0       0
## 5     0       0           0     0 3.188349        1    1          0       0
## 6     0       0           0     0 3.167426        0    0          0       0
men_nofx = subset(men, AnyFx1st == 0)
dim(men_nofx)
## [1] 4909   45

1.2 Men without fracture

1.2.1 Descriptive analysis

createTable(compareGroups(Death ~ age_0 + FNT + B1FND + COMORBIDITIES + cvd + diabetes + copd + cancer + smoke + alcohol_mod + activity, data = men_nofx))
## 
## --------Summary descriptives table by 'Death'---------
## 
## _______________________________________________ 
##                    0           1      p.overall 
##                 N=1991      N=2918              
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## age_0         70.7 (4.62) 75.4 (5.84)  <0.001   
## FNT           0.45 (1.04) 0.59 (1.07)  <0.001   
## B1FND         0.80 (0.13) 0.79 (0.13)  <0.001   
## COMORBIDITIES 0.87 (0.92) 1.37 (1.13)  <0.001   
## cvd           0.13 (0.33) 0.26 (0.44)  <0.001   
## diabetes      0.08 (0.27) 0.13 (0.34)  <0.001   
## copd          0.08 (0.27) 0.12 (0.33)  <0.001   
## cancer        0.14 (0.35) 0.21 (0.41)  <0.001   
## smoke         0.59 (0.49) 0.65 (0.48)  <0.001   
## alcohol_mod   0.25 (0.44) 0.26 (0.44)   0.728   
## activity      0.73 (0.44) 0.63 (0.48)  <0.001   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

1.2.2 Cox’s PH model

Absolute BMD

men_nofx$FNBMD_SD = men_nofx$B1FND/(-0.13)

surv_1 <- Surv(time = men_nofx$entry_end, event = men_nofx$Death)

m.1 <- coxph(surv_1 ~ FNBMD_SD, data = men_nofx)
  summary(m.1)
## Call:
## coxph(formula = surv_1 ~ FNBMD_SD, data = men_nofx)
## 
##   n= 4909, number of events= 2918 
## 
##             coef exp(coef) se(coef)     z Pr(>|z|)    
## FNBMD_SD 0.12275   1.13060  0.01982 6.193 5.91e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## FNBMD_SD     1.131     0.8845     1.088     1.175
## 
## Concordance= 0.539  (se = 0.006 )
## Likelihood ratio test= 39.28  on 1 df,   p=4e-10
## Wald test            = 38.35  on 1 df,   p=6e-10
## Score (logrank) test = 38.3  on 1 df,   p=6e-10
m.1a <- coxph(surv_1 ~ FNBMD_SD + age_0, data = men_nofx)
  summary(m.1a)
## Call:
## coxph(formula = surv_1 ~ FNBMD_SD + age_0, data = men_nofx)
## 
##   n= 4909, number of events= 2918 
## 
##              coef exp(coef) se(coef)      z Pr(>|z|)    
## FNBMD_SD 0.016256  1.016389 0.019365  0.839    0.401    
## age_0    0.121285  1.128947 0.003305 36.698   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## FNBMD_SD     1.016     0.9839    0.9785     1.056
## age_0        1.129     0.8858    1.1217     1.136
## 
## Concordance= 0.685  (se = 0.005 )
## Likelihood ratio test= 1315  on 2 df,   p=<2e-16
## Wald test            = 1395  on 2 df,   p=<2e-16
## Score (logrank) test = 1472  on 2 df,   p=<2e-16
m.1multi <- coxph(surv_1 ~ FNBMD_SD + age_0 + cvd + diabetes + copd + cancer + smoke + alcohol_mod + activity, data = men_nofx)
  summary(m.1multi)
## Call:
## coxph(formula = surv_1 ~ FNBMD_SD + age_0 + cvd + diabetes + 
##     copd + cancer + smoke + alcohol_mod + activity, data = men_nofx)
## 
##   n= 4909, number of events= 2918 
## 
##                  coef exp(coef)  se(coef)      z Pr(>|z|)    
## FNBMD_SD     0.033677  1.034251  0.019483  1.729 0.083889 .  
## age_0        0.116861  1.123964  0.003436 34.014  < 2e-16 ***
## cvd          0.424838  1.529343  0.042890  9.905  < 2e-16 ***
## diabetes     0.459250  1.582887  0.055822  8.227  < 2e-16 ***
## copd         0.270802  1.311016  0.056569  4.787 1.69e-06 ***
## cancer       0.175520  1.191865  0.045838  3.829 0.000129 ***
## smoke        0.292348  1.339569  0.039646  7.374 1.66e-13 ***
## alcohol_mod -0.029633  0.970802  0.042915 -0.691 0.489879    
## activity    -0.190842  0.826263  0.038927 -4.903 9.46e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             exp(coef) exp(-coef) lower .95 upper .95
## FNBMD_SD       1.0343     0.9669    0.9955    1.0745
## age_0          1.1240     0.8897    1.1164    1.1316
## cvd            1.5293     0.6539    1.4060    1.6635
## diabetes       1.5829     0.6318    1.4188    1.7659
## copd           1.3110     0.7628    1.1734    1.4647
## cancer         1.1919     0.8390    1.0895    1.3039
## smoke          1.3396     0.7465    1.2394    1.4478
## alcohol_mod    0.9708     1.0301    0.8925    1.0560
## activity       0.8263     1.2103    0.7656    0.8918
## 
## Concordance= 0.708  (se = 0.005 )
## Likelihood ratio test= 1624  on 9 df,   p=<2e-16
## Wald test            = 1694  on 9 df,   p=<2e-16
## Score (logrank) test = 1804  on 9 df,   p=<2e-16
#ggforest(m.1, data = men_nofx)

BMD T-scores

men_nofx$FNT_SD = men_nofx$FNT/(-1)

m.2 <- coxph(surv_1 ~ FNT_SD, data = men_nofx)
  summary(m.2)
## Call:
## coxph(formula = surv_1 ~ FNT_SD, data = men_nofx)
## 
##   n= 4909, number of events= 2918 
## 
##           coef exp(coef) se(coef)      z Pr(>|z|)    
## FNT_SD -0.1133    0.8929   0.0183 -6.191 5.97e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##        exp(coef) exp(-coef) lower .95 upper .95
## FNT_SD    0.8929       1.12    0.8614    0.9255
## 
## Concordance= 0.539  (se = 0.006 )
## Likelihood ratio test= 39.26  on 1 df,   p=4e-10
## Wald test            = 38.33  on 1 df,   p=6e-10
## Score (logrank) test = 38.28  on 1 df,   p=6e-10
m.2a <- coxph(surv_1 ~ FNT_SD + age_0, data = men_nofx)
  summary(m.2a)
## Call:
## coxph(formula = surv_1 ~ FNT_SD + age_0, data = men_nofx)
## 
##   n= 4909, number of events= 2918 
## 
##             coef exp(coef)  se(coef)      z Pr(>|z|)    
## FNT_SD -0.014956  0.985156  0.017875 -0.837    0.403    
## age_0   0.121287  1.128948  0.003305 36.698   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##        exp(coef) exp(-coef) lower .95 upper .95
## FNT_SD    0.9852     1.0151    0.9512     1.020
## age_0     1.1289     0.8858    1.1217     1.136
## 
## Concordance= 0.685  (se = 0.005 )
## Likelihood ratio test= 1315  on 2 df,   p=<2e-16
## Wald test            = 1395  on 2 df,   p=<2e-16
## Score (logrank) test = 1472  on 2 df,   p=<2e-16
m.2multi <- coxph(surv_1 ~ FNT_SD + age_0 + cvd + diabetes + copd + cancer + smoke + alcohol_mod + activity, data = men_nofx)
  summary(m.2multi)
## Call:
## coxph(formula = surv_1 ~ FNT_SD + age_0 + cvd + diabetes + copd + 
##     cancer + smoke + alcohol_mod + activity, data = men_nofx)
## 
##   n= 4909, number of events= 2918 
## 
##                  coef exp(coef)  se(coef)      z Pr(>|z|)    
## FNT_SD      -0.031051  0.969426  0.017984 -1.727 0.084246 .  
## age_0        0.116862  1.123965  0.003436 34.014  < 2e-16 ***
## cvd          0.424845  1.529353  0.042890  9.905  < 2e-16 ***
## diabetes     0.459237  1.582866  0.055822  8.227  < 2e-16 ***
## copd         0.270804  1.311017  0.056569  4.787 1.69e-06 ***
## cancer       0.175509  1.191852  0.045838  3.829 0.000129 ***
## smoke        0.292350  1.339572  0.039646  7.374 1.66e-13 ***
## alcohol_mod -0.029641  0.970794  0.042915 -0.691 0.489768    
## activity    -0.190846  0.826260  0.038927 -4.903 9.46e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             exp(coef) exp(-coef) lower .95 upper .95
## FNT_SD         0.9694     1.0315    0.9359    1.0042
## age_0          1.1240     0.8897    1.1164    1.1316
## cvd            1.5294     0.6539    1.4060    1.6635
## diabetes       1.5829     0.6318    1.4188    1.7659
## copd           1.3110     0.7628    1.1734    1.4647
## cancer         1.1919     0.8390    1.0894    1.3039
## smoke          1.3396     0.7465    1.2394    1.4478
## alcohol_mod    0.9708     1.0301    0.8925    1.0560
## activity       0.8263     1.2103    0.7656    0.8918
## 
## Concordance= 0.708  (se = 0.005 )
## Likelihood ratio test= 1624  on 9 df,   p=<2e-16
## Wald test            = 1694  on 9 df,   p=<2e-16
## Score (logrank) test = 1804  on 9 df,   p=<2e-16

BMD as a categorical variable

men_nofx$bone[men_nofx$FNT>-1] = "Normal"
men_nofx$bone[men_nofx$FNT>-2.5 & men_nofx$FNT<=-1] = "Osteopenia"
men_nofx$bone[men_nofx$FNT<= -2.5] = "Osteoporosis"

createTable(compareGroups(bone ~ Death + age_0 + FNT + B1FND + COMORBIDITIES + cvd + diabetes + copd + cancer + smoke + alcohol_mod + activity, data = men_nofx))
## 
## --------Summary descriptives table by 'bone'---------
## 
## _____________________________________________________________ 
##                 Normal     Osteopenia  Osteoporosis p.overall 
##                 N=4525       N=349         N=35               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## Death         0.60 (0.49) 0.57 (0.50)  0.60 (0.50)    0.564   
## age_0         73.6 (5.89) 72.0 (5.23)  72.1 (5.45)   <0.001   
## FNT           0.72 (0.85) -1.53 (0.41) -3.37 (0.83)   0.000   
## B1FND         0.77 (0.10) 1.04 (0.05)  1.26 (0.10)    0.000   
## COMORBIDITIES 1.15 (1.08) 1.38 (1.09)  1.37 (1.09)   <0.001   
## cvd           0.20 (0.40) 0.23 (0.42)  0.29 (0.46)    0.213   
## diabetes      0.10 (0.30) 0.19 (0.40)  0.20 (0.41)   <0.001   
## copd          0.11 (0.31) 0.09 (0.28)  0.11 (0.32)    0.589   
## cancer        0.18 (0.39) 0.18 (0.38)  0.20 (0.41)    0.936   
## smoke         0.63 (0.48) 0.62 (0.49)  0.54 (0.51)    0.586   
## alcohol_mod   0.25 (0.43) 0.29 (0.46)  0.29 (0.46)    0.253   
## activity      0.67 (0.47) 0.66 (0.48)  0.69 (0.47)    0.846   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
m.3 <- coxph(surv_1 ~ bone, data = men_nofx)
  summary(m.3)
## Call:
## coxph(formula = surv_1 ~ bone, data = men_nofx)
## 
##   n= 4909, number of events= 2918 
## 
##                      coef exp(coef) se(coef)      z Pr(>|z|)  
## boneOsteopenia   -0.12154   0.88556  0.07363 -1.651   0.0988 .
## boneOsteoporosis -0.08623   0.91738  0.21908 -0.394   0.6939  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                  exp(coef) exp(-coef) lower .95 upper .95
## boneOsteopenia      0.8856      1.129    0.7666     1.023
## boneOsteoporosis    0.9174      1.090    0.5971     1.409
## 
## Concordance= 0.506  (se = 0.002 )
## Likelihood ratio test= 2.95  on 2 df,   p=0.2
## Wald test            = 2.85  on 2 df,   p=0.2
## Score (logrank) test = 2.85  on 2 df,   p=0.2
m.3a <- coxph(surv_1 ~ bone + age_0, data = men_nofx)
summary(m.3a)
## Call:
## coxph(formula = surv_1 ~ bone + age_0, data = men_nofx)
## 
##   n= 4909, number of events= 2918 
## 
##                       coef exp(coef)  se(coef)      z Pr(>|z|)    
## boneOsteopenia    0.026651  1.027009  0.073783  0.361    0.718    
## boneOsteoporosis -0.067972  0.934286  0.219101 -0.310    0.756    
## age_0             0.121812  1.129542  0.003267 37.285   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                  exp(coef) exp(-coef) lower .95 upper .95
## boneOsteopenia      1.0270     0.9737    0.8887     1.187
## boneOsteoporosis    0.9343     1.0703    0.6081     1.435
## age_0               1.1295     0.8853    1.1223     1.137
## 
## Concordance= 0.685  (se = 0.005 )
## Likelihood ratio test= 1315  on 3 df,   p=<2e-16
## Wald test            = 1394  on 3 df,   p=<2e-16
## Score (logrank) test = 1472  on 3 df,   p=<2e-16
m.3multi <- coxph(surv_1 ~ bone + age_0 + cvd + diabetes + copd + cancer + smoke + alcohol_mod + activity, data = men_nofx)
  summary(m.3multi)
## Call:
## coxph(formula = surv_1 ~ bone + age_0 + cvd + diabetes + copd + 
##     cancer + smoke + alcohol_mod + activity, data = men_nofx)
## 
##   n= 4909, number of events= 2918 
## 
##                       coef exp(coef)  se(coef)      z Pr(>|z|)    
## boneOsteopenia   -0.010106  0.989945  0.074150 -0.136 0.891592    
## boneOsteoporosis -0.164170  0.848598  0.219214 -0.749 0.453917    
## age_0             0.117745  1.124957  0.003402 34.610  < 2e-16 ***
## cvd               0.425012  1.529609  0.042889  9.909  < 2e-16 ***
## diabetes          0.447980  1.565148  0.055583  8.060 7.65e-16 ***
## copd              0.273853  1.315022  0.056568  4.841 1.29e-06 ***
## cancer            0.175921  1.192343  0.045837  3.838 0.000124 ***
## smoke             0.292076  1.339205  0.039670  7.363 1.80e-13 ***
## alcohol_mod      -0.032997  0.967541  0.042906 -0.769 0.441851    
## activity         -0.191818  0.825457  0.038924 -4.928 8.31e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                  exp(coef) exp(-coef) lower .95 upper .95
## boneOsteopenia      0.9899     1.0102    0.8560    1.1448
## boneOsteoporosis    0.8486     1.1784    0.5522    1.3041
## age_0               1.1250     0.8889    1.1175    1.1325
## cvd                 1.5296     0.6538    1.4063    1.6638
## diabetes            1.5651     0.6389    1.4036    1.7453
## copd                1.3150     0.7604    1.1770    1.4692
## cancer              1.1923     0.8387    1.0899    1.3044
## smoke               1.3392     0.7467    1.2390    1.4475
## alcohol_mod         0.9675     1.0335    0.8895    1.0524
## activity            0.8255     1.2114    0.7648    0.8909
## 
## Concordance= 0.708  (se = 0.005 )
## Likelihood ratio test= 1621  on 10 df,   p=<2e-16
## Wald test            = 1692  on 10 df,   p=<2e-16
## Score (logrank) test = 1800  on 10 df,   p=<2e-16

1.3 All men

1.3.1 Descriptive analysis

createTable(compareGroups(Death ~ age_0 + FNT + B1FND + COMORBIDITIES + cvd + diabetes + copd + cancer + smoke + alcohol_mod + activity, data = men))
## 
## --------Summary descriptives table by 'Death'---------
## 
## _______________________________________________ 
##                    0           1      p.overall 
##                 N=2382      N=3612              
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## age_0         70.8 (4.66) 75.5 (5.84)  <0.001   
## FNT           0.51 (1.05) 0.68 (1.08)  <0.001   
## B1FND         0.80 (0.13) 0.78 (0.13)  <0.001   
## COMORBIDITIES 0.88 (0.94) 1.36 (1.12)  <0.001   
## cvd           0.13 (0.34) 0.26 (0.44)  <0.001   
## diabetes      0.08 (0.27) 0.13 (0.34)  <0.001   
## copd          0.08 (0.27) 0.13 (0.33)  <0.001   
## cancer        0.14 (0.35) 0.21 (0.41)  <0.001   
## smoke         0.59 (0.49) 0.65 (0.48)  <0.001   
## alcohol_mod   0.25 (0.43) 0.26 (0.44)   0.448   
## activity      0.73 (0.44) 0.62 (0.48)  <0.001   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

1.3.2 Cox’s PH model

Absolute BMD

men$FNBMD_SD = men$B1FND/(-0.13)

surv_2 <- Surv(time = men$entry_end, event = men$Death)

m.21 <- coxph(surv_2 ~ FNBMD_SD, data = men)
  summary(m.21)
## Call:
## coxph(formula = surv_2 ~ FNBMD_SD, data = men)
## 
##   n= 5993, number of events= 3611 
##    (1 observation deleted due to missingness)
## 
##             coef exp(coef) se(coef)     z Pr(>|z|)    
## FNBMD_SD 0.13412   1.14353  0.01778 7.544 4.56e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## FNBMD_SD     1.144     0.8745     1.104     1.184
## 
## Concordance= 0.541  (se = 0.005 )
## Likelihood ratio test= 58.33  on 1 df,   p=2e-14
## Wald test            = 56.91  on 1 df,   p=5e-14
## Score (logrank) test = 56.83  on 1 df,   p=5e-14
m.21a <- coxph(surv_2 ~ FNBMD_SD + age_0, data = men)
  summary(m.21a)
## Call:
## coxph(formula = surv_2 ~ FNBMD_SD + age_0, data = men)
## 
##   n= 5993, number of events= 3611 
##    (1 observation deleted due to missingness)
## 
##              coef exp(coef) se(coef)      z Pr(>|z|)    
## FNBMD_SD 0.020356  1.020564 0.017408  1.169    0.242    
## age_0    0.119908  1.127393 0.002994 40.049   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## FNBMD_SD     1.021     0.9798    0.9863     1.056
## age_0        1.127     0.8870    1.1208     1.134
## 
## 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
m.21multi <- coxph(surv_2 ~ FNBMD_SD + age_0 + cvd + diabetes + copd + cancer + smoke + alcohol_mod + activity, data = men)
  summary(m.21multi)
## Call:
## coxph(formula = surv_2 ~ FNBMD_SD + age_0 + cvd + diabetes + 
##     copd + cancer + smoke + alcohol_mod + activity, data = men)
## 
##   n= 5993, number of events= 3611 
##    (1 observation deleted due to missingness)
## 
##                  coef exp(coef)  se(coef)      z Pr(>|z|)    
## FNBMD_SD     0.037976  1.038706  0.017508  2.169   0.0301 *  
## age_0        0.115895  1.122878  0.003107 37.304  < 2e-16 ***
## cvd          0.411993  1.509823  0.038708 10.644  < 2e-16 ***
## diabetes     0.426617  1.532065  0.050563  8.437  < 2e-16 ***
## copd         0.301771  1.352252  0.050432  5.984 2.18e-09 ***
## cancer       0.174110  1.190186  0.041334  4.212 2.53e-05 ***
## smoke        0.262696  1.300432  0.035588  7.382 1.56e-13 ***
## alcohol_mod -0.016157  0.983973  0.038526 -0.419   0.6749    
## activity    -0.183105  0.832681  0.034975 -5.235 1.65e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             exp(coef) exp(-coef) lower .95 upper .95
## FNBMD_SD       1.0387     0.9627    1.0037    1.0750
## age_0          1.1229     0.8906    1.1161    1.1297
## cvd            1.5098     0.6623    1.3995    1.6288
## diabetes       1.5321     0.6527    1.3875    1.6917
## copd           1.3523     0.7395    1.2250    1.4927
## cancer         1.1902     0.8402    1.0976    1.2906
## smoke          1.3004     0.7690    1.2128    1.3944
## alcohol_mod    0.9840     1.0163    0.9124    1.0611
## activity       0.8327     1.2009    0.7775    0.8918
## 
## Concordance= 0.705  (se = 0.004 )
## Likelihood ratio test= 1938  on 9 df,   p=<2e-16
## Wald test            = 2013  on 9 df,   p=<2e-16
## Score (logrank) test = 2136  on 9 df,   p=<2e-16
#ggforest(m.1, data = men_nofx)

BMD T-scores

men$FNT_SD = men$FNT/(-1)

m.22 <- coxph(surv_2 ~ FNT_SD, data = men)
  summary(m.22)
## Call:
## coxph(formula = surv_2 ~ FNT_SD, data = men)
## 
##   n= 5993, number of events= 3611 
##    (1 observation deleted due to missingness)
## 
##            coef exp(coef) se(coef)      z Pr(>|z|)    
## FNT_SD -0.12376   0.88359  0.01641 -7.541 4.65e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##        exp(coef) exp(-coef) lower .95 upper .95
## FNT_SD    0.8836      1.132    0.8556    0.9125
## 
## Concordance= 0.541  (se = 0.005 )
## Likelihood ratio test= 58.29  on 1 df,   p=2e-14
## Wald test            = 56.87  on 1 df,   p=5e-14
## Score (logrank) test = 56.79  on 1 df,   p=5e-14
m.22a <- coxph(surv_2 ~ FNT_SD + age_0, data = men)
  summary(m.22a)
## Call:
## coxph(formula = surv_2 ~ FNT_SD + age_0, data = men)
## 
##   n= 5993, number of events= 3611 
##    (1 observation deleted due to missingness)
## 
##             coef exp(coef)  se(coef)      z Pr(>|z|)    
## FNT_SD -0.018731  0.981443  0.016069 -1.166    0.244    
## age_0   0.119909  1.127395  0.002994 40.049   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##        exp(coef) exp(-coef) lower .95 upper .95
## FNT_SD    0.9814      1.019     0.951     1.013
## age_0     1.1274      0.887     1.121     1.134
## 
## 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
m.22multi <- coxph(surv_2 ~ FNT_SD + age_0 + cvd + diabetes + copd + cancer + smoke + alcohol_mod + activity, data = men)
  summary(m.22multi)
## Call:
## coxph(formula = surv_2 ~ FNT_SD + age_0 + cvd + diabetes + copd + 
##     cancer + smoke + alcohol_mod + activity, data = men)
## 
##   n= 5993, number of events= 3611 
##    (1 observation deleted due to missingness)
## 
##                  coef exp(coef)  se(coef)      z Pr(>|z|)    
## FNT_SD      -0.035002  0.965604  0.016161 -2.166   0.0303 *  
## age_0        0.115896  1.122880  0.003107 37.304  < 2e-16 ***
## cvd          0.411997  1.509830  0.038708 10.644  < 2e-16 ***
## diabetes     0.426595  1.532033  0.050563  8.437  < 2e-16 ***
## copd         0.301772  1.352253  0.050432  5.984 2.18e-09 ***
## cancer       0.174104  1.190180  0.041334  4.212 2.53e-05 ***
## smoke        0.262698  1.300434  0.035588  7.382 1.56e-13 ***
## alcohol_mod -0.016168  0.983962  0.038526 -0.420   0.6747    
## activity    -0.183107  0.832679  0.034975 -5.235 1.65e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             exp(coef) exp(-coef) lower .95 upper .95
## FNT_SD         0.9656     1.0356    0.9355    0.9967
## age_0          1.1229     0.8906    1.1161    1.1297
## cvd            1.5098     0.6623    1.3995    1.6288
## diabetes       1.5320     0.6527    1.3875    1.6916
## copd           1.3523     0.7395    1.2250    1.4927
## cancer         1.1902     0.8402    1.0976    1.2906
## smoke          1.3004     0.7690    1.2128    1.3944
## alcohol_mod    0.9840     1.0163    0.9124    1.0611
## activity       0.8327     1.2009    0.7775    0.8918
## 
## Concordance= 0.705  (se = 0.004 )
## Likelihood ratio test= 1938  on 9 df,   p=<2e-16
## Wald test            = 2013  on 9 df,   p=<2e-16
## Score (logrank) test = 2136  on 9 df,   p=<2e-16

BMD as a categorical variable

men$bone[men$FNT>-1] = "Normal"
men$bone[men$FNT>-2.5 & men$FNT<=-1] = "Osteopenia"
men$bone[men$FNT<= -2.5] = "Osteoporosis"

createTable(compareGroups(bone ~ Death + age_0 + FNT + B1FND + COMORBIDITIES + cvd + diabetes + copd + cancer + smoke + alcohol_mod + activity, data = men))
## 
## --------Summary descriptives table by 'bone'---------
## 
## _____________________________________________________________ 
##                 Normal     Osteopenia  Osteoporosis p.overall 
##                 N=5570       N=387         N=36               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## Death         0.61 (0.49) 0.57 (0.50)  0.58 (0.50)    0.355   
## age_0         73.8 (5.90) 72.0 (5.19)  72.2 (5.39)   <0.001   
## FNT           0.79 (0.87) -1.53 (0.41) -3.37 (0.82)   0.000   
## B1FND         0.76 (0.10) 1.04 (0.05)  1.26 (0.10)    0.000   
## COMORBIDITIES 1.16 (1.08) 1.36 (1.09)  1.39 (1.08)    0.001   
## cvd           0.21 (0.41) 0.23 (0.42)  0.31 (0.47)    0.225   
## diabetes      0.10 (0.30) 0.19 (0.39)  0.19 (0.40)   <0.001   
## copd          0.11 (0.31) 0.09 (0.29)  0.11 (0.32)    0.662   
## cancer        0.18 (0.39) 0.17 (0.37)  0.19 (0.40)    0.665   
## smoke         0.62 (0.48) 0.64 (0.48)  0.56 (0.50)    0.627   
## alcohol_mod   0.25 (0.44) 0.31 (0.46)  0.28 (0.45)    0.054   
## activity      0.67 (0.47) 0.65 (0.48)  0.69 (0.47)    0.708   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
m.23 <- coxph(surv_2 ~ bone, data = men)
  summary(m.23)
## Call:
## coxph(formula = surv_2 ~ bone, data = men)
## 
##   n= 5993, number of events= 3611 
##    (1 observation deleted due to missingness)
## 
##                      coef exp(coef) se(coef)      z Pr(>|z|)  
## boneOsteopenia   -0.12487   0.88261  0.06959 -1.794   0.0728 .
## boneOsteoporosis -0.11851   0.88824  0.21891 -0.541   0.5882  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                  exp(coef) exp(-coef) lower .95 upper .95
## boneOsteopenia      0.8826      1.133    0.7701     1.012
## boneOsteoporosis    0.8882      1.126    0.5784     1.364
## 
## Concordance= 0.505  (se = 0.002 )
## Likelihood ratio test= 3.6  on 2 df,   p=0.2
## Wald test            = 3.48  on 2 df,   p=0.2
## Score (logrank) test = 3.48  on 2 df,   p=0.2
m.23a <- coxph(surv_2 ~ bone + age_0, data = men)
summary(m.23a)
## Call:
## coxph(formula = surv_2 ~ bone + age_0, data = men)
## 
##   n= 5993, number of events= 3611 
##    (1 observation deleted due to missingness)
## 
##                       coef exp(coef)  se(coef)      z Pr(>|z|)    
## boneOsteopenia    0.046450  1.047546  0.069762  0.666    0.506    
## boneOsteoporosis -0.079803  0.923298  0.218921 -0.365    0.715    
## age_0             0.120658  1.128239  0.002954 40.847   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                  exp(coef) exp(-coef) lower .95 upper .95
## boneOsteopenia      1.0475     0.9546    0.9137     1.201
## boneOsteoporosis    0.9233     1.0831    0.6012     1.418
## age_0               1.1282     0.8863    1.1217     1.135
## 
## Concordance= 0.682  (se = 0.005 )
## Likelihood ratio test= 1585  on 3 df,   p=<2e-16
## Wald test            = 1673  on 3 df,   p=<2e-16
## Score (logrank) test = 1761  on 3 df,   p=<2e-16
m.23multi <- coxph(surv_2 ~ bone + age_0 + cvd + diabetes + copd + cancer + smoke + alcohol_mod + activity, data = men)
  summary(m.23multi)
## Call:
## coxph(formula = surv_2 ~ bone + age_0 + cvd + diabetes + copd + 
##     cancer + smoke + alcohol_mod + activity, data = men)
## 
##   n= 5993, number of events= 3611 
##    (1 observation deleted due to missingness)
## 
##                       coef exp(coef)  se(coef)      z Pr(>|z|)    
## boneOsteopenia    0.013728  1.013823  0.070040  0.196    0.845    
## boneOsteoporosis -0.181585  0.833948  0.219044 -0.829    0.407    
## age_0             0.117028  1.124151  0.003073 38.086  < 2e-16 ***
## cvd               0.411826  1.509571  0.038714 10.638  < 2e-16 ***
## diabetes          0.413146  1.511565  0.050340  8.207 2.27e-16 ***
## copd              0.305714  1.357593  0.050412  6.064 1.33e-09 ***
## cancer            0.174641  1.190818  0.041337  4.225 2.39e-05 ***
## smoke             0.262249  1.299850  0.035599  7.367 1.75e-13 ***
## alcohol_mod      -0.022318  0.977929  0.038492 -0.580    0.562    
## activity         -0.183165  0.832631  0.034976 -5.237 1.63e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                  exp(coef) exp(-coef) lower .95 upper .95
## boneOsteopenia      1.0138     0.9864    0.8838    1.1630
## boneOsteoporosis    0.8339     1.1991    0.5429    1.2811
## age_0               1.1242     0.8896    1.1174    1.1309
## cvd                 1.5096     0.6624    1.3993    1.6286
## diabetes            1.5116     0.6616    1.3695    1.6683
## copd                1.3576     0.7366    1.2299    1.4986
## cancer              1.1908     0.8398    1.0981    1.2913
## smoke               1.2998     0.7693    1.2122    1.3938
## alcohol_mod         0.9779     1.0226    0.9069    1.0546
## activity            0.8326     1.2010    0.7775    0.8917
## 
## Concordance= 0.705  (se = 0.004 )
## Likelihood ratio test= 1934  on 10 df,   p=<2e-16
## Wald test            = 2009  on 10 df,   p=<2e-16
## Score (logrank) test = 2131  on 10 df,   p=<2e-16

2. WOMEN (SOF data)

2.1 Data setup

women = read.csv("C:\\UTS\\Research projects\\Skeletal Age_BMD and Mortality\\SOF data.csv")
head(women)
##   ID FX50 intercept        slope       Rate Base.FNBMD DEATH AGE WEIGHT HEIGHT
## 1  1    1 0.6696049 -0.002579665 -0.3932416      0.656     0  69   67.3  150.5
## 2  3    0 0.6533052 -0.008523688 -1.3033162      0.654     0  75   72.7  151.0
## 3  4    0 0.8238722  0.001288358  0.1530117      0.842     0  67   60.6  155.1
## 4  5    0 0.4977875 -0.011105498 -2.2803898      0.487     0  84   44.7  157.0
## 5 10    1 0.6035814  0.004313417  0.7249440      0.595     0  69   85.2  165.1
## 6 12    1 0.7023159 -0.002350381 -0.3451366      0.681     1  75   71.0  156.4
##        BMI SMOKE Drink Timetodeath        BMIG Smoke  Ageg Drinking Death
## 1 29.71270     1     1    17.23562  Overweight  Past 64-69      Yes     0
## 2 31.88457     1     1    18.64932       Obese  Past 75-79      Yes     0
## 3 25.19121     0     1    19.76164  Overweight Never 64-69      Yes     0
## 4 18.13461     0     1    14.88219 Underweight Never   80+      Yes     0
## 5 31.25687     0     0    17.23562       Obese Never 64-69       No     0
## 6 29.02584     0     1    19.36164  Overweight Never 75-79      Yes     1
##   Quartile_Rate Group_boneloss   AGE_SD FNBMD_SD   BMI_SD     Rate_1
## 1             3          Least 15.82569 6.074074 7.125347  0.3932416
## 2             1           Most 17.20183 6.055556 7.646179  1.3033162
## 3             4          Least 15.36697 7.796296 6.041058 -0.1530117
## 4             1           Most 19.26606 4.509259 4.348827  2.2803898
## 5             4          Least 15.82569 5.509259 7.495652 -0.7249440
## 6             3          Least 17.20183 6.305556 6.960634  0.3451366
##          RateG EST BISH CAL
## 1      -1 to 1  No   No Yes
## 2     -2 to -1  No  Yes Yes
## 3      -1 to 1  No  Yes Yes
## 4 less than -2  No   No Yes
## 5      -1 to 1  No   No Yes
## 6      -1 to 1  No   No Yes
#women_nofx = subset(women, AnyFx1st == 0)
#dim(women_nofx)

2.2 All women

2.2.1 Descriptive analysis

createTable(compareGroups(Death ~ AGE + Base.FNBMD + SMOKE + Drink, data = women))
## 
## --------Summary descriptives table by 'Death'---------
## 
## ____________________________________________ 
##                 0           1      p.overall 
##              N=1881      N=1942              
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## AGE        71.0 (3.55) 73.6 (4.69)  <0.001   
## Base.FNBMD 0.66 (0.11) 0.65 (0.11)  <0.001   
## SMOKE      0.48 (0.62) 0.57 (0.68)  <0.001   
## Drink      0.81 (0.39) 0.78 (0.41)   0.036   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

1.2.2 Cox’s PH model

Absolute BMD

women$FNBMD_SD = women$Base.FNBMD/(-0.11)

surv_1w <- Surv(time = women$Timetodeath, event = women$Death)

w.1 <- coxph(surv_1w ~ FNBMD_SD, data = women)
  summary(w.1)
## Call:
## coxph(formula = surv_1w ~ FNBMD_SD, data = women)
## 
##   n= 3823, number of events= 1942 
## 
##             coef exp(coef) se(coef)     z Pr(>|z|)    
## FNBMD_SD 0.13807   1.14805  0.02439 5.662  1.5e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## FNBMD_SD     1.148      0.871     1.094     1.204
## 
## Concordance= 0.545  (se = 0.007 )
## Likelihood ratio test= 33.09  on 1 df,   p=9e-09
## Wald test            = 32.06  on 1 df,   p=1e-08
## Score (logrank) test = 31.99  on 1 df,   p=2e-08
w.1a <- coxph(surv_1w ~ FNBMD_SD + AGE, data = women)
  summary(w.1a)
## Call:
## coxph(formula = surv_1w ~ FNBMD_SD + AGE, data = women)
## 
##   n= 3823, number of events= 1942 
## 
##              coef exp(coef) se(coef)      z Pr(>|z|)    
## FNBMD_SD 0.016170  1.016301 0.024432  0.662    0.508    
## AGE      0.124361  1.132425 0.005065 24.555   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## FNBMD_SD     1.016     0.9840    0.9688     1.066
## AGE          1.132     0.8831    1.1212     1.144
## 
## Concordance= 0.648  (se = 0.007 )
## Likelihood ratio test= 571.2  on 2 df,   p=<2e-16
## Wald test            = 652.3  on 2 df,   p=<2e-16
## Score (logrank) test = 679  on 2 df,   p=<2e-16
w.1multi <- coxph(surv_1w ~ FNBMD_SD + AGE + SMOKE + Drink, data = women)
  summary(w.1multi)
## Call:
## coxph(formula = surv_1w ~ FNBMD_SD + AGE + SMOKE + Drink, data = women)
## 
##   n= 3823, number of events= 1942 
## 
##              coef exp(coef) se(coef)      z Pr(>|z|)    
## FNBMD_SD  0.01167   1.01174  0.02438  0.479   0.6321    
## AGE       0.12785   1.13639  0.00512 24.974  < 2e-16 ***
## SMOKE     0.24570   1.27851  0.03447  7.128 1.02e-12 ***
## Drink    -0.12499   0.88251  0.05509 -2.269   0.0233 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## FNBMD_SD    1.0117     0.9884    0.9645    1.0613
## AGE         1.1364     0.8800    1.1250    1.1478
## SMOKE       1.2785     0.7822    1.1950    1.3679
## Drink       0.8825     1.1331    0.7922    0.9831
## 
## Concordance= 0.655  (se = 0.006 )
## Likelihood ratio test= 622.7  on 4 df,   p=<2e-16
## Wald test            = 695.8  on 4 df,   p=<2e-16
## Score (logrank) test = 725.9  on 4 df,   p=<2e-16
#ggforest(m.1, data = men_nofx)