Association between BMD and mortality risk: A multistate approach

Analysis plan

Rationales: - Previous studies have shown a significantly positive association betwen BMD and mortality risk. Nevertheless, the occurance of incident fractures was not accounted for in the previous studies. Given a proven mortality risk following fragility fractures, it is critical to control for the occurence of incident fractures in any assessment of the association between BMD and mortality risk. - Evidence of “off-label” treatment of osteoporosis on mortality has been found, though there is no accepted biomechanisms. - Robust evidence of beneficial contribution of BMD (i.e., ~ possible benefits of osteoporosis treatments) would widen the treatment indications, ultimately reducing the risk of mortality. Analysis approach: - Multistate regression to differentiate the contribution of BMD to (i) incident fracture (which is well established), (ii) death without fractures, and (iii) post-fracture deaths.

Data management

bmd.0 = read.csv("C:\\Thach\\Research projects\\BMD and mortality\\Data\\SOF_MrOS_Nick_21mar24.csv")

bmd = subset(bmd.0, select = c("ID", "gender", "age", "race", "education", "weight", "height", "BMI", "smoke", "drink", "fall", "fx50", "physical", "hypertension", "copd", "parkinson", "cancer", "rheumatoid", "cvd", "renal", "depression", "diabetes", "fnbmd", "anyfx", "death", "Tscore", "event", "state", "time2event", "ageBase", "fnbmdBase", "TscoreBase", "time2end"))
dim(bmd)
## [1] 30538    33
head(bmd)
##   ID gender   age    race education weight height     BMI smoke drink fall fx50
## 1  1      F 69.00 1:WHITE         9   67.3  150.5 29.7127     1     1    2    1
## 2  1      F 86.01 1:WHITE         9   59.5  148.8 26.8727     1     1    2    1
## 3  2      F 84.00 1:WHITE        10   68.4  151.8 29.6833     0     0    1    1
## 4  2      F 90.53 1:WHITE        10   68.4  151.8 29.6833     0     0    1    1
## 5  3      F 75.00 1:WHITE        12   72.7  151.0 31.8846     1     1    3    0
## 6  3      F 93.46 1:WHITE        12   61.7  147.2 28.4754     1     1    3    0
##   physical hypertension copd parkinson cancer rheumatoid cvd renal depression
## 1        0            0    0         0      0          0   0     0          0
## 2        0            0    0         0      0          0   0     0          0
## 3        0            1    0         0      0          0   1     0          0
## 4        0            1    0         0      0          0   1     0          0
## 5        1            0    1         0      1          0   0     0          0
## 6        1            0    1         0      1          0   0     0          0
##   diabetes  fnbmd anyfx death  Tscore event state time2event ageBase fnbmdBase
## 1        0 0.6560     0     0 -1.7000  Well     1     0.0000      69     0.656
## 2        0 0.6349     0     0 -1.8758  Well     1    17.2238      69     0.656
## 3        2 0.5850     0     1 -2.2917  Well     1     0.0000      84     0.585
## 4        2 0.5850     0     1 -2.2917 Death     3     6.5325      84     0.585
## 5        0 0.6540     0     0 -1.7167  Well     1     0.0000      75     0.654
## 6        0 0.5537     0     0 -2.5525  Well     1    18.6366      75     0.654
##   TscoreBase time2end
## 1    -1.7000  17.2238
## 2    -1.7000   5.0103
## 3    -2.2917   6.5325
## 4    -2.2917   6.5325
## 5    -1.7167  18.6366
## 6    -1.7167   6.4586
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.3     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
bmd = bmd %>% mutate(fx_death = case_when(anyfx == 0 & death == 0 ~ "Event-free",
                                          anyfx == 1 & death == 0 ~ "Fractured, Alive",
                                          anyfx == 0 & death == 1 ~ "No Fractured, Dead",
                                          anyfx == 1 & death == 1 ~ "Fractured, Dead"))
bmd = bmd %>% mutate(fall.no = case_when(fall == 0 ~ "0",
                                         fall == 1 ~ "1",
                                         fall == 2 ~ "2",
                                         fall >= 3 ~ "3+"))
bmd = bmd %>% mutate(fall.yesno = case_when(fall == 0 ~ "No",
                                            fall >= 1 ~ "Yes"))
bmd$Tscore.2 = bmd$Tscore*(-1)
bmd$TscoreBase.2 = bmd$TscoreBase*(-1)
bmd = bmd %>% mutate(cvd.n = case_when(cvd == 0 ~ "No",
                                       cvd >= 1 ~ "Yes"))
bmd = bmd %>% mutate(diabetes.n = case_when(diabetes == 0 ~ "No",
                                            diabetes >= 1 ~ "Yes"))
bmd = bmd %>% mutate(drink.n = case_when(drink == 0 ~ "No",
                                         drink >= 1 ~ "Yes"))

1. MEN

1.1 Description of study sample

men = subset(bmd, gender == "M" & race == "1:WHITE")
baseline.m = subset(men, state == 1 & time2event == 0)
library(table1)
## 
## Attaching package: 'table1'
## The following objects are masked from 'package:base':
## 
##     units, units<-
table1(~ ageBase + fnbmdBase + TscoreBase + as.factor(fall.no) + fall.yesno + as.factor(fx50) + race + weight + 
         height + BMI + as.factor(smoke) + as.factor(drink.n) + as.factor(physical) + as.factor(cvd.n) +
         as.factor(hypertension) + as.factor(copd) + as.factor(diabetes.n) + as.factor(cancer) + as.factor(parkinson) +
         as.factor(rheumatoid) + as.factor(renal) + as.factor(depression) + as.factor(anyfx) + as.factor(death) |
         fx_death, data = baseline.m)
Event-free
(N=1722)
Fractured, Alive
(N=361)
Fractured, Dead
(N=655)
No Fractured, Dead
(N=2641)
Overall
(N=5379)
ageBase
Mean (SD) 70.8 (4.66) 71.7 (4.87) 76.2 (5.83) 75.6 (5.87) 73.8 (5.92)
Median [Min, Max] 70.0 [64.0, 89.0] 71.0 [65.0, 91.0] 76.0 [65.0, 94.0] 75.0 [65.0, 100] 73.0 [64.0, 100]
fnbmdBase
Mean (SD) 0.803 (0.122) 0.754 (0.115) 0.727 (0.121) 0.782 (0.122) 0.780 (0.124)
Median [Min, Max] 0.794 [0.404, 1.32] 0.744 [0.499, 1.27] 0.719 [0.273, 1.15] 0.773 [0.475, 1.33] 0.772 [0.273, 1.33]
TscoreBase
Mean (SD) -0.929 (0.894) -1.29 (0.843) -1.48 (0.883) -1.08 (0.893) -1.09 (0.905)
Median [Min, Max] -0.992 [-3.84, 2.84] -1.36 [-3.15, 2.49] -1.54 [-4.80, 1.59] -1.15 [-3.32, 2.89] -1.16 [-4.80, 2.89]
as.factor(fall.no)
0 1344 (78.0%) 248 (68.7%) 389 (59.4%) 1966 (74.4%) 3947 (73.4%)
1 378 (22.0%) 113 (31.3%) 266 (40.6%) 675 (25.6%) 1432 (26.6%)
fall.yesno
No 1344 (78.0%) 248 (68.7%) 389 (59.4%) 1966 (74.4%) 3947 (73.4%)
Yes 378 (22.0%) 113 (31.3%) 266 (40.6%) 675 (25.6%) 1432 (26.6%)
as.factor(fx50)
0 1475 (85.7%) 280 (77.6%) 492 (75.1%) 2186 (82.8%) 4433 (82.4%)
1 247 (14.3%) 81 (22.4%) 163 (24.9%) 455 (17.2%) 946 (17.6%)
race
1:WHITE 1722 (100%) 361 (100%) 655 (100%) 2641 (100%) 5379 (100%)
weight
Mean (SD) 84.4 (12.2) 83.9 (12.7) 82.3 (13.3) 83.2 (13.5) 83.5 (13.0)
Median [Min, Max] 83.3 [52.7, 129] 82.6 [58.0, 142] 80.8 [52.6, 142] 81.5 [50.8, 144] 82.1 [50.8, 144]
height
Mean (SD) 175 (6.41) 176 (6.90) 174 (6.69) 174 (6.66) 174 (6.63)
Median [Min, Max] 175 [154, 198] 175 [158, 199] 174 [153, 193] 174 [147, 198] 174 [147, 199]
Missing 1 (0.1%) 1 (0.3%) 0 (0%) 4 (0.2%) 6 (0.1%)
BMI
Mean (SD) 27.5 (3.56) 27.2 (3.56) 27.1 (3.87) 27.4 (3.96) 27.4 (3.80)
Median [Min, Max] 27.1 [17.5, 50.7] 26.6 [19.2, 41.9] 26.6 [18.3, 41.5] 26.8 [17.2, 48.5] 26.9 [17.2, 50.7]
Missing 1 (0.1%) 1 (0.3%) 0 (0%) 4 (0.2%) 6 (0.1%)
as.factor(smoke)
0 708 (41.1%) 143 (39.6%) 245 (37.4%) 918 (34.8%) 2014 (37.4%)
1 1014 (58.9%) 218 (60.4%) 410 (62.6%) 1723 (65.2%) 3365 (62.6%)
as.factor(drink.n)
No 547 (31.8%) 111 (30.7%) 238 (36.3%) 945 (35.8%) 1841 (34.2%)
Yes 1175 (68.2%) 250 (69.3%) 417 (63.7%) 1696 (64.2%) 3538 (65.8%)
as.factor(physical)
0 458 (26.6%) 96 (26.6%) 255 (38.9%) 976 (37.0%) 1785 (33.2%)
1 1264 (73.4%) 265 (73.4%) 400 (61.1%) 1665 (63.0%) 3594 (66.8%)
as.factor(cvd.n)
No 1507 (87.5%) 304 (84.2%) 497 (75.9%) 1946 (73.7%) 4254 (79.1%)
Yes 215 (12.5%) 57 (15.8%) 158 (24.1%) 695 (26.3%) 1125 (20.9%)
as.factor(hypertension)
0 1113 (64.6%) 231 (64.0%) 344 (52.5%) 1436 (54.4%) 3124 (58.1%)
1 609 (35.4%) 130 (36.0%) 311 (47.5%) 1205 (45.6%) 2255 (41.9%)
as.factor(copd)
0 1583 (91.9%) 334 (92.5%) 562 (85.8%) 2321 (87.9%) 4800 (89.2%)
1 139 (8.1%) 27 (7.5%) 93 (14.2%) 320 (12.1%) 579 (10.8%)
as.factor(diabetes.n)
No 1611 (93.6%) 338 (93.6%) 582 (88.9%) 2313 (87.6%) 4844 (90.1%)
Yes 111 (6.4%) 23 (6.4%) 73 (11.1%) 328 (12.4%) 535 (9.9%)
as.factor(cancer)
0 1481 (86.0%) 305 (84.5%) 529 (80.8%) 2077 (78.6%) 4392 (81.7%)
1 241 (14.0%) 56 (15.5%) 126 (19.2%) 564 (21.4%) 987 (18.3%)
as.factor(parkinson)
0 1505 (87.4%) 311 (86.1%) 579 (88.4%) 2254 (85.3%) 4649 (86.4%)
1 217 (12.6%) 50 (13.9%) 76 (11.6%) 387 (14.7%) 730 (13.6%)
as.factor(rheumatoid)
0 982 (57.0%) 201 (55.7%) 335 (51.1%) 1314 (49.8%) 2832 (52.6%)
1 740 (43.0%) 160 (44.3%) 320 (48.9%) 1327 (50.2%) 2547 (47.4%)
as.factor(renal)
0 1512 (87.8%) 307 (85.0%) 599 (91.5%) 2446 (92.6%) 4864 (90.4%)
1 210 (12.2%) 54 (15.0%) 56 (8.5%) 195 (7.4%) 515 (9.6%)
as.factor(depression)
0 1658 (96.3%) 345 (95.6%) 622 (95.0%) 2551 (96.6%) 5176 (96.2%)
1 64 (3.7%) 16 (4.4%) 33 (5.0%) 90 (3.4%) 203 (3.8%)
as.factor(anyfx)
0 1722 (100%) 0 (0%) 0 (0%) 2641 (100%) 4363 (81.1%)
1 0 (0%) 361 (100%) 655 (100%) 0 (0%) 1016 (18.9%)
as.factor(death)
0 1722 (100%) 361 (100%) 0 (0%) 0 (0%) 2083 (38.7%)
1 0 (0%) 0 (0%) 655 (100%) 2641 (100%) 3296 (61.3%)
#library(compareGroups)
#createTable(compareGroups(fx_death ~ ageBase + fnbmdBase + TscoreBase + fall.no + fall.yesno + fx50 + weight + height + BMI + smoke + drink.n + physical + cvd.n + hypertension + copd + diabetes.n + cancer + parkinson + rheumatoid + renal + depression + anyfx + death, data = baseline.m))

1.2 Conventional Cox’s PH model

1.2.1 Age-adjusted model
library(survival)
cox.m = coxph(Surv(time2end, death) ~ ageBase + TscoreBase.2, data = baseline.m)
summary(cox.m)
## Call:
## coxph(formula = Surv(time2end, death) ~ ageBase + TscoreBase.2, 
##     data = baseline.m)
## 
##   n= 5379, number of events= 3296 
## 
##                  coef exp(coef) se(coef)      z Pr(>|z|)    
## ageBase      0.119424  1.126848 0.003122 38.252   <2e-16 ***
## TscoreBase.2 0.043428  1.044385 0.019810  2.192   0.0284 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##              exp(coef) exp(-coef) lower .95 upper .95
## ageBase          1.127     0.8874     1.120     1.134
## TscoreBase.2     1.044     0.9575     1.005     1.086
## 
## Concordance= 0.683  (se = 0.005 )
## Likelihood ratio test= 1465  on 2 df,   p=<2e-16
## Wald test            = 1544  on 2 df,   p=<2e-16
## Score (logrank) test = 1624  on 2 df,   p=<2e-16
1.2.3 Adjusted for essential covariates (~ JBMR paper)
cox.m2 = coxph(Surv(time2end, death) ~ ageBase + TscoreBase.2 + fall.yesno + fx50 + cvd.n + copd + diabetes.n + cancer, data = baseline.m)
summary(cox.m2)
## Call:
## coxph(formula = Surv(time2end, death) ~ ageBase + TscoreBase.2 + 
##     fall.yesno + fx50 + cvd.n + copd + diabetes.n + cancer, data = baseline.m)
## 
##   n= 5379, number of events= 3296 
## 
##                    coef exp(coef)  se(coef)      z Pr(>|z|)    
## ageBase        0.116148  1.123162  0.003213 36.149  < 2e-16 ***
## TscoreBase.2   0.056079  1.057682  0.019974  2.808 0.004991 ** 
## fall.yesnoYes -0.051258  0.950034  0.039096 -1.311 0.189829    
## fx50           0.061496  1.063426  0.045091  1.364 0.172620    
## cvd.nYes       0.411689  1.509365  0.040635 10.131  < 2e-16 ***
## copd           0.344438  1.411196  0.052879  6.514 7.33e-11 ***
## diabetes.nYes  0.481182  1.617985  0.054185  8.880  < 2e-16 ***
## cancer         0.154837  1.167468  0.043331  3.573 0.000352 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## ageBase           1.123     0.8903    1.1161     1.130
## TscoreBase.2      1.058     0.9455    1.0171     1.100
## fall.yesnoYes     0.950     1.0526    0.8800     1.026
## fx50              1.063     0.9404    0.9735     1.162
## cvd.nYes          1.509     0.6625    1.3938     1.634
## copd              1.411     0.7086    1.2723     1.565
## diabetes.nYes     1.618     0.6181    1.4550     1.799
## cancer            1.167     0.8566    1.0724     1.271
## 
## Concordance= 0.7  (se = 0.005 )
## Likelihood ratio test= 1705  on 8 df,   p=<2e-16
## Wald test            = 1776  on 8 df,   p=<2e-16
## Score (logrank) test = 1884  on 8 df,   p=<2e-16
1.2.3 Multivariable-adjusted model
cox.m3 = coxph(Surv(time2end, death) ~ ageBase + TscoreBase.2 + fall.yesno + fx50 + BMI + smoke + drink.n + 
                 physical + cvd.n + hypertension + copd + diabetes.n + cancer + renal + parkinson + 
                 depression, data = baseline.m)
summary(cox.m3)
## Call:
## coxph(formula = Surv(time2end, death) ~ ageBase + TscoreBase.2 + 
##     fall.yesno + fx50 + BMI + smoke + drink.n + physical + cvd.n + 
##     hypertension + copd + diabetes.n + cancer + renal + parkinson + 
##     depression, data = baseline.m)
## 
##   n= 5373, number of events= 3292 
##    (6 observations deleted due to missingness)
## 
##                    coef exp(coef)  se(coef)       z Pr(>|z|)    
## ageBase        0.114047  1.120805  0.003318  34.375  < 2e-16 ***
## TscoreBase.2   0.075446  1.078365  0.021064   3.582 0.000341 ***
## fall.yesnoYes -0.052975  0.948403  0.039254  -1.350 0.177154    
## fx50           0.067608  1.069946  0.045136   1.498 0.134166    
## BMI            0.013421  1.013512  0.005303   2.531 0.011375 *  
## smoke          0.243538  1.275755  0.037240   6.540 6.16e-11 ***
## drink.nYes    -0.097462  0.907137  0.037266  -2.615 0.008914 ** 
## physical      -0.162434  0.850072  0.036899  -4.402 1.07e-05 ***
## cvd.nYes       0.385894  1.470928  0.040914   9.432  < 2e-16 ***
## hypertension   0.180788  1.198161  0.036079   5.011 5.42e-07 ***
## copd           0.276813  1.318919  0.053335   5.190 2.10e-07 ***
## diabetes.nYes  0.397446  1.488020  0.055429   7.170 7.48e-13 ***
## cancer         0.131224  1.140224  0.043562   3.012 0.002592 ** 
## renal         -0.878005  0.415611  0.080458 -10.913  < 2e-16 ***
## parkinson      0.532339  1.702911  0.061545   8.650  < 2e-16 ***
## depression    -0.112735  0.893387  0.092733  -1.216 0.224100    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## ageBase          1.1208     0.8922    1.1135    1.1281
## TscoreBase.2     1.0784     0.9273    1.0348    1.1238
## fall.yesnoYes    0.9484     1.0544    0.8782    1.0242
## fx50             1.0699     0.9346    0.9794    1.1689
## BMI              1.0135     0.9867    1.0030    1.0241
## smoke            1.2758     0.7838    1.1860    1.3724
## drink.nYes       0.9071     1.1024    0.8432    0.9759
## physical         0.8501     1.1764    0.7908    0.9138
## cvd.nYes         1.4709     0.6798    1.3576    1.5937
## hypertension     1.1982     0.8346    1.1164    1.2860
## copd             1.3189     0.7582    1.1880    1.4643
## diabetes.nYes    1.4880     0.6720    1.3348    1.6588
## cancer           1.1402     0.8770    1.0469    1.2419
## renal            0.4156     2.4061    0.3550    0.4866
## parkinson        1.7029     0.5872    1.5094    1.9212
## depression       0.8934     1.1193    0.7449    1.0715
## 
## Concordance= 0.718  (se = 0.004 )
## Likelihood ratio test= 1956  on 16 df,   p=<2e-16
## Wald test            = 2017  on 16 df,   p=<2e-16
## Score (logrank) test = 2132  on 16 df,   p=<2e-16

1.3 Multistate regression

1.3.1 Data preparation

reorder_array <- function(Qx3) {
  new_row_order <- c("Well", "Fracture", "Death")
  new_col_order <- c("Well", "Fracture")
  Qx3_reordered <- Qx3[new_col_order, new_row_order]
  Death <- c(0, 0, 0)
  Qx3_reordered <- rbind(Qx3_reordered, Death)
  return(Qx3_reordered)
}
library(msm)
## Warning: package 'msm' was built under R version 4.3.3
Q.m = statetable.msm(event, subject = ID, data = bmd)
Q.m2 = reorder_array(Q.m)
qmatrix.n = crudeinits.msm(state ~ time2event, subject = ID, data = men, qmatrix = Q.m2)
qmatrix.n
##                 Well   Fracture      Death
## Well     -0.05500799  0.0152825 0.03972549
## Fracture  0.00000000 -0.1033188 0.10331884
## Death     0.00000000  0.0000000 0.00000000
statetable.msm(state, ID, data = men)
##     to
## from    1    2    3
##    1 1722 1016 2641
##    2    0  361  655

1.3.2 Age-adjusted model:

Age and FNBMD at baseline:
age.m1= msm(state ~ time2event, subject = ID, data = men, qmatrix = qmatrix.n, gen.inits = TRUE, exacttimes = TRUE, 
           method = "BFGS", control = list(fnscale = 4000, maxit = 10000), pci = 5, covariates =~ ageBase + TscoreBase.2)  
age.m1
## 
## Call:
## msm(formula = state ~ time2event, subject = ID, data = men, qmatrix = qmatrix.n,     gen.inits = TRUE, covariates = ~ageBase + TscoreBase.2, exacttimes = TRUE,     pci = 5, method = "BFGS", control = list(fnscale = 4000,         maxit = 10000))
## 
## Maximum likelihood estimates
## Baselines are with covariates set to their means
## 
## Transition intensities with hazard ratios for each covariate
##                     Baseline                     ageBase            
## Well - Well         -0.04267 (-0.04452,-0.04089)                    
## Well - Fracture      0.01365 ( 0.01271, 0.01467) 1.061 (1.049,1.072)
## Well - Death         0.02901 ( 0.02753, 0.03058) 1.117 (1.109,1.124)
## Fracture - Fracture -0.07334 (-0.08547,-0.06294)                    
## Fracture - Death     0.07334 ( 0.06294, 0.08547) 1.103 (1.088,1.119)
##                     TscoreBase.2          timeperiod[5,Inf)  
## Well - Well                                                  
## Well - Fracture     1.6307 (1.5095,1.762) 1.541 (1.347,1.763)
## Well - Death        0.9809 (0.9399,1.024) 3.308 (2.997,3.651)
## Fracture - Fracture                                          
## Fracture - Death    1.1391 (1.0396,1.248) 1.657 (1.258,2.182)
## 
## -2 * log-likelihood:  35082.2
hazard.msm(age.m1)
## $ageBase
##                        HR        L        U
## Well - Fracture  1.060684 1.049403 1.072086
## Well - Death     1.116725 1.109467 1.124031
## Fracture - Death 1.103430 1.088470 1.118597
## 
## $TscoreBase.2
##                        HR         L        U
## Well - Fracture  1.630696 1.5095166 1.761603
## Well - Death     0.980934 0.9398526 1.023811
## Fracture - Death 1.139145 1.0395698 1.248258
## 
## $`timeperiod[5,Inf)`
##                        HR        L        U
## Well - Fracture  1.541144 1.347192 1.763018
## Well - Death     3.307866 2.996597 3.651467
## Fracture - Death 1.657043 1.258378 2.182008
Age and FNBMD at fracture time:
age.m2= msm(state ~ time2event, subject = ID, data = men, qmatrix = qmatrix.n, gen.inits = TRUE, exacttimes = TRUE, 
           method = "BFGS", control = list(fnscale = 4000, maxit = 10000), pci = 5, covariates =~ age + Tscore.2)  
age.m2
## 
## Call:
## msm(formula = state ~ time2event, subject = ID, data = men, qmatrix = qmatrix.n,     gen.inits = TRUE, covariates = ~age + Tscore.2, exacttimes = TRUE,     pci = 5, method = "BFGS", control = list(fnscale = 4000,         maxit = 10000))
## 
## Maximum likelihood estimates
## Baselines are with covariates set to their means
## 
## Transition intensities with hazard ratios for each covariate
##                     Baseline                     age                
## Well - Well         -0.04604 (-0.04799,-0.04417)                    
## Well - Fracture      0.01432 ( 0.01334, 0.01538) 1.061 (1.049,1.072)
## Well - Death         0.03172 ( 0.03014, 0.03338) 1.117 (1.109,1.124)
## Fracture - Fracture -0.05318 (-0.06310,-0.04481)                    
## Fracture - Death     0.05318 ( 0.04481, 0.06310) 1.094 (1.081,1.108)
##                     Tscore.2              timeperiod[5,Inf)    
## Well - Well                                                    
## Well - Fracture     1.6308 (1.5096,1.762) 1.5407 (1.3468,1.762)
## Well - Death        0.9809 (0.9399,1.024) 3.3078 (2.9965,3.651)
## Fracture - Fracture                                            
## Fracture - Death    1.1555 (1.0536,1.267) 0.9588 (0.7299,1.259)
## 
## -2 * log-likelihood:  35069.19
hazard.msm(age.m2)
## $age
##                        HR        L        U
## Well - Fracture  1.060684 1.049404 1.072086
## Well - Death     1.116699 1.109441 1.124005
## Fracture - Death 1.094149 1.080765 1.107698
## 
## $Tscore.2
##                        HR         L        U
## Well - Fracture  1.630793 1.5096138 1.761699
## Well - Death     0.980941 0.9398584 1.023819
## Fracture - Death 1.155495 1.0535774 1.267271
## 
## $`timeperiod[5,Inf)`
##                         HR         L        U
## Well - Fracture  1.5406553 1.3467891 1.762428
## Well - Death     3.3077907 2.9965212 3.651394
## Fracture - Death 0.9587954 0.7298911 1.259487

1.3.3 Adjusted for essential covariates:

Age and FNBMD at baseline:
multi.me1= msm(state ~ time2event, subject = ID, data = men, qmatrix = qmatrix.n, gen.inits = TRUE, exacttimes = TRUE, 
           method = "BFGS", control = list(fnscale = 4000, maxit = 10000), pci = 5, covariates =~ ageBase + TscoreBase.2 +
             fall.yesno + fx50 + cvd.n + hypertension + copd + diabetes.n + cancer)  
multi.me1
## 
## Call:
## msm(formula = state ~ time2event, subject = ID, data = men, qmatrix = qmatrix.n,     gen.inits = TRUE, covariates = ~ageBase + TscoreBase.2 +         fall.yesno + fx50 + cvd.n + hypertension + copd + diabetes.n +         cancer, exacttimes = TRUE, pci = 5, method = "BFGS",     control = list(fnscale = 4000, maxit = 10000))
## 
## Maximum likelihood estimates
## Baselines are with covariates set to their means
## 
## Transition intensities with hazard ratios for each covariate
##                     Baseline                     ageBase            
## Well - Well         -0.04168 (-0.04352,-0.03992)                    
## Well - Fracture      0.01337 ( 0.01243, 0.01438) 1.055 (1.043,1.067)
## Well - Death         0.02831 ( 0.02684, 0.02986) 1.112 (1.104,1.119)
## Fracture - Fracture -0.07164 (-0.08373,-0.06129)                    
## Fracture - Death     0.07164 ( 0.06129, 0.08373) 1.102 (1.087,1.118)
##                     TscoreBase.2         fall.yesnoYes         
## Well - Well                                                    
## Well - Fracture     1.634 (1.5104,1.767) 1.6040 (1.4107,1.8237)
## Well - Death        1.001 (0.9587,1.045) 0.8754 (0.8014,0.9563)
## Fracture - Fracture                                            
## Fracture - Death    1.153 (1.0513,1.265) 0.9155 (0.7799,1.0747)
##                     fx50                  cvd.nYes            
## Well - Well                                                   
## Well - Fracture     1.3744 (1.1882,1.590) 1.145 (0.9817,1.336)
## Well - Death        1.0308 (0.9308,1.141) 1.512 (1.3839,1.652)
## Fracture - Fracture                                           
## Fracture - Death    0.9982 (0.8336,1.195) 1.173 (0.9758,1.409)
##                     hypertension        copd                
## Well - Well                                                 
## Well - Fracture     1.198 (1.056,1.358) 1.143 (0.9435,1.384)
## Well - Death        1.187 (1.098,1.283) 1.333 (1.1849,1.499)
## Fracture - Fracture                                         
## Fracture - Death    1.319 (1.127,1.545) 1.478 (1.1844,1.845)
##                     diabetes.nYes       cancer               
## Well - Well                                                  
## Well - Fracture     1.272 (1.028,1.575) 0.9571 (0.8135,1.126)
## Well - Death        1.564 (1.390,1.760) 1.1817 (1.0753,1.299)
## Fracture - Fracture                                          
## Fracture - Death    1.293 (1.003,1.668) 1.0148 (0.8333,1.236)
##                     timeperiod[5,Inf)  
## Well - Well                            
## Well - Fracture     1.594 (1.393,1.824)
## Well - Death        3.458 (3.132,3.819)
## Fracture - Fracture                    
## Fracture - Death    1.724 (1.307,2.274)
## 
## -2 * log-likelihood:  34736.31
hazard.msm(multi.me1)
## $ageBase
##                        HR        L        U
## Well - Fracture  1.054914 1.043404 1.066552
## Well - Death     1.111779 1.104302 1.119306
## Fracture - Death 1.102462 1.087109 1.118031
## 
## $TscoreBase.2
##                        HR         L        U
## Well - Fracture  1.633550 1.5103854 1.766758
## Well - Death     1.000990 0.9586538 1.045195
## Fracture - Death 1.153403 1.0513061 1.265414
## 
## $fall.yesnoYes
##                         HR         L         U
## Well - Fracture  1.6039588 1.4106632 1.8237406
## Well - Death     0.8754452 0.8014228 0.9563046
## Fracture - Death 0.9154675 0.7798582 1.0746579
## 
## $fx50
##                         HR         L        U
## Well - Fracture  1.3743877 1.1882493 1.589684
## Well - Death     1.0307597 0.9308102 1.141442
## Fracture - Death 0.9982193 0.8336076 1.195337
## 
## $cvd.nYes
##                        HR         L        U
## Well - Fracture  1.145045 0.9816527 1.335633
## Well - Death     1.512191 1.3838808 1.652399
## Fracture - Death 1.172537 0.9758125 1.408922
## 
## $hypertension
##                        HR        L        U
## Well - Fracture  1.197531 1.055667 1.358458
## Well - Death     1.186891 1.098088 1.282876
## Fracture - Death 1.319439 1.126843 1.544953
## 
## $copd
##                        HR        L        U
## Well - Fracture  1.142676 0.943541 1.383839
## Well - Death     1.332596 1.184945 1.498645
## Fracture - Death 1.478103 1.184388 1.844655
## 
## $diabetes.nYes
##                        HR        L        U
## Well - Fracture  1.272375 1.028001 1.574842
## Well - Death     1.563967 1.389978 1.759736
## Fracture - Death 1.293216 1.002616 1.668045
## 
## $cancer
##                        HR         L        U
## Well - Fracture  0.957080 0.8135196 1.125974
## Well - Death     1.181678 1.0752816 1.298603
## Fracture - Death 1.014791 0.8333236 1.235775
## 
## $`timeperiod[5,Inf)`
##                        HR        L        U
## Well - Fracture  1.593833 1.392623 1.824114
## Well - Death     3.458342 3.132120 3.818541
## Fracture - Death 1.723676 1.306746 2.273632
Age and FNBMD at fracture time:
multi.me2= msm(state ~ time2event, subject = ID, data = men, qmatrix = qmatrix.n, gen.inits = TRUE, exacttimes = TRUE, 
           method = "BFGS", control = list(fnscale = 4000, maxit = 10000), pci = 5, covariates =~ age + Tscore.2 +
             fall.yesno + fx50 + cvd.n + hypertension + copd + diabetes.n + cancer)  
multi.me2
## 
## Call:
## msm(formula = state ~ time2event, subject = ID, data = men, qmatrix = qmatrix.n,     gen.inits = TRUE, covariates = ~age + Tscore.2 + fall.yesno +         fx50 + cvd.n + hypertension + copd + diabetes.n + cancer,     exacttimes = TRUE, pci = 5, method = "BFGS", control = list(fnscale = 4000,         maxit = 10000))
## 
## Maximum likelihood estimates
## Baselines are with covariates set to their means
## 
## Transition intensities with hazard ratios for each covariate
##                     Baseline                     age                
## Well - Well         -0.04480 (-0.04674,-0.04295)                    
## Well - Fracture      0.01396 ( 0.01299, 0.01501) 1.055 (1.043,1.067)
## Well - Death         0.03084 ( 0.02928, 0.03249) 1.112 (1.104,1.119)
## Fracture - Fracture -0.04902 (-0.05858,-0.04102)                    
## Fracture - Death     0.04902 ( 0.04102, 0.05858) 1.097 (1.083,1.110)
##                     Tscore.2             fall.yesnoYes         
## Well - Well                                                    
## Well - Fracture     1.634 (1.5106,1.767) 1.6032 (1.4100,1.8229)
## Well - Death        1.001 (0.9586,1.045) 0.8754 (0.8014,0.9562)
## Fracture - Fracture                                            
## Fracture - Death    1.172 (1.0671,1.286) 1.0560 (0.8995,1.2397)
##                     fx50                 cvd.nYes            
## Well - Well                                                  
## Well - Fracture     1.374 (1.1883,1.590) 1.145 (0.9818,1.336)
## Well - Death        1.031 (0.9308,1.141) 1.512 (1.3838,1.652)
## Fracture - Fracture                                          
## Fracture - Death    1.097 (0.9148,1.315) 1.266 (1.0546,1.520)
##                     hypertension        copd                
## Well - Well                                                 
## Well - Fracture     1.197 (1.056,1.358) 1.143 (0.9438,1.384)
## Well - Death        1.187 (1.098,1.283) 1.333 (1.1852,1.499)
## Fracture - Fracture                                         
## Fracture - Death    1.315 (1.122,1.541) 1.578 (1.2637,1.969)
##                     diabetes.nYes       cancer               
## Well - Well                                                  
## Well - Fracture     1.272 (1.028,1.575) 0.9569 (0.8133,1.126)
## Well - Death        1.564 (1.390,1.760) 1.1818 (1.0754,1.299)
## Fracture - Fracture                                          
## Fracture - Death    1.361 (1.055,1.756) 1.0482 (0.8611,1.276)
##                     timeperiod[5,Inf)   
## Well - Well                             
## Well - Fracture     1.594 (1.3924,1.824)
## Well - Death        3.458 (3.1318,3.818)
## Fracture - Fracture                     
## Fracture - Death    1.035 (0.7859,1.362)
## 
## -2 * log-likelihood:  34709.41
hazard.msm(multi.me2)
## $age
##                        HR        L        U
## Well - Fracture  1.054908 1.043398 1.066545
## Well - Death     1.111766 1.104289 1.119293
## Fracture - Death 1.096603 1.082895 1.110484
## 
## $Tscore.2
##                        HR         L        U
## Well - Fracture  1.633751 1.5105733 1.766973
## Well - Death     1.000972 0.9586371 1.045177
## Fracture - Death 1.171577 1.0671066 1.286275
## 
## $fall.yesnoYes
##                         HR         L         U
## Well - Fracture  1.6032429 1.4100305 1.8229307
## Well - Death     0.8753802 0.8013617 0.9562354
## Fracture - Death 1.0560090 0.8995297 1.2397090
## 
## $fx50
##                        HR         L        U
## Well - Fracture  1.374487 1.1883459 1.589786
## Well - Death     1.030766 0.9308155 1.141449
## Fracture - Death 1.096847 0.9148146 1.315100
## 
## $cvd.nYes
##                        HR         L        U
## Well - Fracture  1.145187 0.9817847 1.335786
## Well - Death     1.512109 1.3838021 1.652313
## Fracture - Death 1.266245 1.0546011 1.520362
## 
## $hypertension
##                        HR        L        U
## Well - Fracture  1.197473 1.055621 1.358386
## Well - Death     1.186594 1.097812 1.282556
## Fracture - Death 1.314628 1.121581 1.540903
## 
## $copd
##                        HR        L        U
## Well - Fracture  1.142981 0.943816 1.384173
## Well - Death     1.332856 1.185185 1.498927
## Fracture - Death 1.577621 1.263716 1.969500
## 
## $diabetes.nYes
##                        HR        L        U
## Well - Fracture  1.272274 1.027918 1.574718
## Well - Death     1.564041 1.390043 1.759819
## Fracture - Death 1.360782 1.054662 1.755753
## 
## $cancer
##                         HR         L        U
## Well - Fracture  0.9568559 0.8133231 1.125719
## Well - Death     1.1817624 1.0753587 1.298694
## Fracture - Death 1.0482357 0.8610708 1.276083
## 
## $`timeperiod[5,Inf)`
##                        HR         L        U
## Well - Fracture  1.593514 1.3923577 1.823732
## Well - Death     3.457999 3.1318135 3.818157
## Fracture - Death 1.034763 0.7858665 1.362488

1.3.4 Multivariable-adjusted model:

All covariates at baseline:
multi.m1 = msm(state ~ time2event, subject = ID, data = men, qmatrix = qmatrix.n, gen.inits = TRUE, exacttimes = TRUE, 
           method = "BFGS", control = list(fnscale = 4000, maxit = 10000), pci = 5, covariates =~ ageBase + TscoreBase.2 +
             fall.yesno + fx50 + BMI + smoke + drink.n + physical + cvd.n + hypertension + copd + diabetes.n + cancer + 
             renal + parkinson + depression)  
multi.m1
## 
## Call:
## msm(formula = state ~ time2event, subject = ID, data = men, qmatrix = qmatrix.n,     gen.inits = TRUE, covariates = ~ageBase + TscoreBase.2 +         fall.yesno + fx50 + BMI + smoke + drink.n + physical +         cvd.n + hypertension + copd + diabetes.n + cancer + renal +         parkinson + depression, exacttimes = TRUE, pci = 5, method = "BFGS",     control = list(fnscale = 4000, maxit = 10000))
## 
## Maximum likelihood estimates
## Baselines are with covariates set to their means
## 
## Transition intensities with hazard ratios for each covariate
##                     Baseline                     ageBase            
## Well - Well         -0.04108 (-0.04291,-0.03933)                    
## Well - Fracture      0.01334 ( 0.01240, 0.01436) 1.057 (1.045,1.069)
## Well - Death         0.02774 ( 0.02628, 0.02928) 1.108 (1.100,1.116)
## Fracture - Fracture -0.06969 (-0.08165,-0.05949)                    
## Fracture - Death     0.06969 ( 0.05949, 0.08165) 1.105 (1.089,1.121)
##                     TscoreBase.2         fall.yesnoYes         
## Well - Well                                                    
## Well - Fracture     1.675 (1.5428,1.819) 1.5993 (1.4059,1.8193)
## Well - Death        1.010 (0.9653,1.057) 0.8769 (0.8025,0.9582)
## Fracture - Fracture                                            
## Fracture - Death    1.183 (1.0726,1.304) 0.8717 (0.7407,1.0259)
##                     fx50                  BMI                
## Well - Well                                                  
## Well - Fracture     1.3650 (1.1797,1.579) 1.020 (1.002,1.039)
## Well - Death        1.0261 (0.9265,1.136) 1.013 (1.001,1.025)
## Fracture - Fracture                                          
## Fracture - Death    0.9874 (0.8242,1.183) 1.007 (0.984,1.031)
##                     smoke                drink.nYes            
## Well - Well                                                    
## Well - Fracture     1.067 (0.9381,1.214) 0.9768 (0.8559,1.1148)
## Well - Death        1.275 (1.1749,1.383) 0.9084 (0.8373,0.9856)
## Fracture - Fracture                                            
## Fracture - Death    1.183 (1.0034,1.394) 0.9248 (0.7826,1.0929)
##                     physical               cvd.nYes            
## Well - Well                                                    
## Well - Fracture     0.9272 (0.8124,1.0581) 1.133 (0.9703,1.323)
## Well - Death        0.8606 (0.7937,0.9331) 1.492 (1.3653,1.632)
## Fracture - Fracture                                            
## Fracture - Death    0.9165 (0.7776,1.0803) 1.156 (0.9614,1.389)
##                     hypertension        copd                
## Well - Well                                                 
## Well - Fracture     1.166 (1.026,1.325) 1.123 (0.9258,1.361)
## Well - Death        1.142 (1.056,1.236) 1.254 (1.1140,1.412)
## Fracture - Fracture                                         
## Fracture - Death    1.309 (1.115,1.538) 1.488 (1.1885,1.864)
##                     diabetes.nYes        cancer               
## Well - Well                                                   
## Well - Fracture     1.232 (0.9933,1.529) 0.9546 (0.8111,1.123)
## Well - Death        1.492 (1.3232,1.681) 1.1795 (1.0730,1.297)
## Fracture - Fracture                                           
## Fracture - Death    1.281 (0.9880,1.661) 0.9774 (0.8004,1.194)
##                     renal                  parkinson            
## Well - Well                                                     
## Well - Fracture     1.0167 (0.7851,1.3166) 0.8601 (0.6750,1.096)
## Well - Death        0.4328 (0.3624,0.5168) 1.7169 (1.5057,1.958)
## Fracture - Fracture                                             
## Fracture - Death    0.4589 (0.3315,0.6354) 1.1814 (0.8887,1.570)
##                     depression            timeperiod[5,Inf)  
## Well - Well                                                  
## Well - Fracture     1.0617 (0.7943,1.419) 1.606 (1.403,1.839)
## Well - Death        0.8303 (0.6718,1.026) 3.621 (3.277,3.999)
## Fracture - Fracture                                          
## Fracture - Death    1.2174 (0.8500,1.743) 1.857 (1.404,2.457)
## 
## -2 * log-likelihood:  34482.87
hazard.msm(multi.m1)
## $ageBase
##                        HR        L        U
## Well - Fracture  1.056546 1.044640 1.068587
## Well - Death     1.108130 1.100438 1.115877
## Fracture - Death 1.104521 1.088564 1.120712
## 
## $TscoreBase.2
##                        HR         L        U
## Well - Fracture  1.675293 1.5428062 1.819157
## Well - Death     1.010163 0.9653202 1.057088
## Fracture - Death 1.182816 1.0726030 1.304353
## 
## $fall.yesnoYes
##                         HR         L         U
## Well - Fracture  1.5993353 1.4059333 1.8193420
## Well - Death     0.8769080 0.8024858 0.9582321
## Fracture - Death 0.8717051 0.7406923 1.0258913
## 
## $fx50
##                         HR         L        U
## Well - Fracture  1.3650267 1.1797272 1.579431
## Well - Death     1.0260644 0.9265252 1.136297
## Fracture - Death 0.9873546 0.8241698 1.182850
## 
## $BMI
##                        HR         L        U
## Well - Fracture  1.020471 1.0018024 1.039488
## Well - Death     1.012820 1.0012579 1.024515
## Fracture - Death 1.006990 0.9839527 1.030567
## 
## $smoke
##                        HR         L        U
## Well - Fracture  1.067370 0.9381202 1.214427
## Well - Death     1.274774 1.1748640 1.383181
## Fracture - Death 1.182677 1.0034061 1.393977
## 
## $drink.nYes
##                         HR         L         U
## Well - Fracture  0.9768408 0.8559284 1.1148338
## Well - Death     0.9084156 0.8372520 0.9856279
## Fracture - Death 0.9248342 0.7825822 1.0929436
## 
## $physical
##                         HR         L         U
## Well - Fracture  0.9271842 0.8124486 1.0581231
## Well - Death     0.8605990 0.7937252 0.9331071
## Fracture - Death 0.9165490 0.7776290 1.0802866
## 
## $cvd.nYes
##                        HR         L        U
## Well - Fracture  1.132837 0.9703376 1.322551
## Well - Death     1.492496 1.3653119 1.631527
## Fracture - Death 1.155582 0.9614215 1.388953
## 
## $hypertension
##                        HR        L        U
## Well - Fracture  1.166079 1.026074 1.325187
## Well - Death     1.142402 1.055782 1.236128
## Fracture - Death 1.309404 1.115075 1.537599
## 
## $copd
##                        HR        L        U
## Well - Fracture  1.122663 0.925805 1.361380
## Well - Death     1.254072 1.113987 1.411771
## Fracture - Death 1.488416 1.188492 1.864027
## 
## $diabetes.nYes
##                        HR         L        U
## Well - Fracture  1.232288 0.9932675 1.528827
## Well - Death     1.491572 1.3231849 1.681387
## Fracture - Death 1.280897 0.9879516 1.660706
## 
## $cancer
##                         HR         L        U
## Well - Fracture  0.9545850 0.8110881 1.123469
## Well - Death     1.1795137 1.0730453 1.296546
## Fracture - Death 0.9774277 0.8004478 1.193538
## 
## $renal
##                         HR         L         U
## Well - Fracture  1.0167070 0.7851332 1.3165832
## Well - Death     0.4328100 0.3624450 0.5168356
## Fracture - Death 0.4589377 0.3314948 0.6353758
## 
## $parkinson
##                         HR         L        U
## Well - Fracture  0.8601457 0.6750162 1.096049
## Well - Death     1.7169301 1.5056776 1.957822
## Fracture - Death 1.1813560 0.8886750 1.570430
## 
## $depression
##                        HR         L        U
## Well - Fracture  1.061699 0.7943326 1.419058
## Well - Death     0.830307 0.6717878 1.026232
## Fracture - Death 1.217374 0.8500413 1.743444
## 
## $`timeperiod[5,Inf)`
##                        HR        L        U
## Well - Fracture  1.606265 1.402872 1.839145
## Well - Death     3.620526 3.277497 3.999457
## Fracture - Death 1.857024 1.403578 2.456963
Age and FNBMD at fracture time (other covariates at baseline):
multi.m2= msm(state ~ time2event, subject = ID, data = men, qmatrix = qmatrix.n, gen.inits = TRUE, exacttimes = TRUE, 
           method = "BFGS", control = list(fnscale = 4000, maxit = 10000), pci = 5, covariates =~ age + Tscore.2 + fall.yesno
           + fx50 + BMI + smoke + drink + physical + cvd + hypertension + copd + diabetes + cancer + renal + parkinson 
           + depression)  
multi.m2
## 
## Call:
## msm(formula = state ~ time2event, subject = ID, data = men, qmatrix = qmatrix.n,     gen.inits = TRUE, covariates = ~age + Tscore.2 + fall.yesno +         fx50 + BMI + smoke + drink + physical + cvd + hypertension +         copd + diabetes + cancer + renal + parkinson + depression,     exacttimes = TRUE, pci = 5, method = "BFGS", control = list(fnscale = 4000,         maxit = 10000))
## 
## Maximum likelihood estimates
## Baselines are with covariates set to their means
## 
## Transition intensities with hazard ratios for each covariate
##                     Baseline                     age                
## Well - Well         -0.04410 (-0.04603,-0.04226)                    
## Well - Fracture      0.01395 ( 0.01297, 0.01500) 1.057 (1.045,1.069)
## Well - Death         0.03015 ( 0.02860, 0.03178) 1.108 (1.101,1.116)
## Fracture - Fracture -0.04795 (-0.05740,-0.04006)                    
## Fracture - Death     0.04795 ( 0.04006, 0.05740) 1.096 (1.082,1.110)
##                     Tscore.2             fall.yesnoYes         
## Well - Well                                                    
## Well - Fracture     1.675 (1.5421,1.819) 1.6002 (1.4065,1.8204)
## Well - Death        1.012 (0.9671,1.059) 0.8747 (0.8004,0.9558)
## Fracture - Fracture                                            
## Fracture - Death    1.205 (1.0903,1.331) 1.0177 (0.8653,1.1969)
##                     fx50                 BMI                 
## Well - Well                                                  
## Well - Fracture     1.365 (1.1796,1.579) 1.020 (1.0017,1.039)
## Well - Death        1.027 (0.9271,1.137) 1.013 (1.0018,1.025)
## Fracture - Fracture                                          
## Fracture - Death    1.085 (0.9047,1.302) 1.010 (0.9861,1.034)
##                     smoke                drink                
## Well - Well                                                   
## Well - Fracture     1.069 (0.9383,1.217) 0.9930 (0.9509,1.037)
## Well - Death        1.268 (1.1677,1.376) 0.9889 (0.9626,1.016)
## Fracture - Fracture                                           
## Fracture - Death    1.135 (0.9626,1.339) 0.9953 (0.9406,1.053)
##                     physical               cvd                 
## Well - Well                                                    
## Well - Fracture     0.9271 (0.8124,1.0581) 1.133 (0.9706,1.323)
## Well - Death        0.8594 (0.7926,0.9318) 1.500 (1.3720,1.639)
## Fracture - Fracture                                            
## Fracture - Death    0.8557 (0.7277,1.0063) 1.247 (1.0370,1.499)
##                     hypertension        copd                
## Well - Well                                                 
## Well - Fracture     1.166 (1.026,1.326) 1.123 (0.9262,1.362)
## Well - Death        1.143 (1.056,1.237) 1.260 (1.1192,1.418)
## Fracture - Fracture                                         
## Fracture - Death    1.302 (1.109,1.530) 1.552 (1.2396,1.943)
##                     diabetes             cancer               
## Well - Well                                                   
## Well - Fracture     1.233 (0.9937,1.529) 0.9542 (0.8108,1.123)
## Well - Death        1.505 (1.3351,1.696) 1.1770 (1.0708,1.294)
## Fracture - Fracture                                           
## Fracture - Death    1.341 (1.0342,1.739) 1.0073 (0.8252,1.230)
##                     renal                  parkinson            
## Well - Well                                                     
## Well - Fracture     1.0173 (0.7856,1.3173) 0.8598 (0.6747,1.096)
## Well - Death        0.4316 (0.3614,0.5154) 1.7245 (1.5122,1.967)
## Fracture - Fracture                                             
## Fracture - Death    0.4689 (0.3351,0.6560) 1.2100 (0.9017,1.624)
##                     depression            timeperiod[5,Inf)   
## Well - Well                                                   
## Well - Fracture     1.0600 (0.7932,1.417) 1.606 (1.4026,1.839)
## Well - Death        0.8277 (0.6697,1.023) 3.615 (3.2727,3.994)
## Fracture - Fracture                                           
## Fracture - Death    1.2079 (0.8453,1.726) 1.112 (0.8412,1.469)
## 
## -2 * log-likelihood:  34462.53
hazard.msm(multi.m2)
## $age
##                        HR        L        U
## Well - Fracture  1.056570 1.044664 1.068611
## Well - Death     1.108328 1.100629 1.116082
## Fracture - Death 1.096260 1.082355 1.110344
## 
## $Tscore.2
##                        HR         L        U
## Well - Fracture  1.674721 1.5420523 1.818803
## Well - Death     1.012175 0.9671222 1.059327
## Fracture - Death 1.204762 1.0902598 1.331290
## 
## $fall.yesnoYes
##                         HR         L         U
## Well - Fracture  1.6001564 1.4065287 1.8204396
## Well - Death     0.8746621 0.8004058 0.9558075
## Fracture - Death 1.0176745 0.8652604 1.1969360
## 
## $fx50
##                        HR         L        U
## Well - Fracture  1.364864 1.1795886 1.579240
## Well - Death     1.026662 0.9270577 1.136967
## Fracture - Death 1.085207 0.9047267 1.301691
## 
## $BMI
##                        HR         L        U
## Well - Fracture  1.020409 1.0017204 1.039447
## Well - Death     1.013402 1.0018080 1.025130
## Fracture - Death 1.009760 0.9860529 1.034036
## 
## $smoke
##                        HR         L        U
## Well - Fracture  1.068511 0.9382613 1.216841
## Well - Death     1.267724 1.1677324 1.376278
## Fracture - Death 1.135320 0.9626011 1.339030
## 
## $drink
##                         HR         L        U
## Well - Fracture  0.9930048 0.9508621 1.037015
## Well - Death     0.9888865 0.9626160 1.015874
## Fracture - Death 0.9952812 0.9405593 1.053187
## 
## $physical
##                         HR         L        U
## Well - Fracture  0.9271293 0.8124043 1.058055
## Well - Death     0.8594246 0.7926377 0.931839
## Fracture - Death 0.8557376 0.7277214 1.006274
## 
## $cvd
##                        HR         L        U
## Well - Fracture  1.133114 0.9706489 1.322771
## Well - Death     1.499757 1.3719823 1.639432
## Fracture - Death 1.246841 1.0370190 1.499116
## 
## $hypertension
##                        HR        L        U
## Well - Fracture  1.166422 1.026394 1.325555
## Well - Death     1.142921 1.056257 1.236694
## Fracture - Death 1.302444 1.108841 1.529849
## 
## $copd
##                        HR         L        U
## Well - Fracture  1.123038 0.9261638 1.361760
## Well - Death     1.259936 1.1192194 1.418344
## Fracture - Death 1.551957 1.2396025 1.943019
## 
## $diabetes
##                        HR        L        U
## Well - Fracture  1.232738 0.993742 1.529214
## Well - Death     1.504881 1.335069 1.696291
## Fracture - Death 1.341176 1.034180 1.739304
## 
## $cancer
##                         HR         L        U
## Well - Fracture  0.9541689 0.8107589 1.122946
## Well - Death     1.1770365 1.0707994 1.293814
## Fracture - Death 1.0072977 0.8251633 1.229634
## 
## $renal
##                         HR         L         U
## Well - Fracture  1.0172769 0.7855852 1.3173012
## Well - Death     0.4315937 0.3614036 0.5154156
## Fracture - Death 0.4688591 0.3351264 0.6559582
## 
## $parkinson
##                        HR         L        U
## Well - Fracture  0.859820 0.6747487 1.095653
## Well - Death     1.724504 1.5121765 1.966645
## Fracture - Death 1.210036 0.9017124 1.623785
## 
## $depression
##                         HR         L        U
## Well - Fracture  1.0599914 0.7931759 1.416561
## Well - Death     0.8277142 0.6696892 1.023028
## Fracture - Death 1.2078658 0.8453157 1.725911
## 
## $`timeperiod[5,Inf)`
##                        HR         L        U
## Well - Fracture  1.605877 1.4025580 1.838669
## Well - Death     3.615157 3.2726543 3.993504
## Fracture - Death 1.111629 0.8411597 1.469066

1.3.4 Transition risks

# Osteopenia:
options(digits = 5)
pmatrix.msm(multi.m2, t = 1, ci = "none", covariates = list(age = 60, Tscore.2 = -1.5), 2)
##             Well   Fracture     Death
## Well     0.99744 0.00056862 0.0019935
## Fracture 0.00000 0.99524036 0.0047596
## Death    0.00000 0.00000000 1.0000000
  pmatrix.msm(multi.m2, t = 5, ci = "none", covariates = list(age = 60, Tscore.2 = -1.5), 2)
##             Well  Fracture    Death
## Well     0.97633 0.0034702 0.020197
## Fracture 0.00000 0.9753878 0.024612
## Death    0.00000 0.0000000 1.000000
  pmatrix.msm(multi.m2, t = 10, ci = "none", covariates = list(age = 60, Tscore.2 = -1.5), 2)
##             Well  Fracture    Death
## Well     0.93745 0.0077057 0.054843
## Fracture 0.00000 0.9498625 0.050138
## Death    0.00000 0.0000000 1.000000
# Osteoporosis:
pmatrix.msm(multi.m2, t = 1, ci = "none", covariates = list(age = 60, Tscore.2 = -2.5), 2)
##             Well   Fracture     Death
## Well     0.99769 0.00033971 0.0019691
## Fracture 0.00000 0.99604771 0.0039523
## Death    0.00000 0.00000000 1.0000000
  pmatrix.msm(multi.m2, t = 5, ci = "none", covariates = list(age = 60, Tscore.2 = -2.5), 2)
##             Well  Fracture    Death
## Well     0.97797 0.0020778 0.019954
## Fracture 0.00000 0.9795277 0.020472
## Death    0.00000 0.0000000 1.000000
  pmatrix.msm(multi.m2, t = 10, ci = "none", covariates = list(age = 60, Tscore.2 = -2.5), 2)
##             Well Fracture    Death
## Well     0.94117 0.004629 0.054206
## Fracture 0.00000 0.958203 0.041797
## Death    0.00000 0.000000 1.000000
# Osteoporosis + comorbidities:
pmatrix.msm(multi.m2, t = 1, ci = "none", covariates = list(age = 60, Tscore.2 = -2.5, cvd = 1, diabetes = 1, cancer = 1), 2)
##             Well  Fracture     Death
## Well     0.99433 0.0004514 0.0052219
## Fracture 0.00000 0.9933516 0.0066484
## Death    0.00000 0.0000000 1.0000000
  pmatrix.msm(multi.m2, t = 5, ci = "none", covariates = list(age = 60, Tscore.2 = -2.5, cvd = 1, diabetes = 1, cancer = 1), 2)
##             Well  Fracture    Death
## Well     0.94518 0.0027137 0.052106
## Fracture 0.00000 0.9657579 0.034242
## Death    0.00000 0.0000000 1.000000
  pmatrix.msm(multi.m2, t = 10, ci = "none", covariates = list(age = 60, Tscore.2 = -2.5, cvd = 1, diabetes = 1, cancer = 1), 2)
##             Well  Fracture    Death
## Well     0.85671 0.0058369 0.137455
## Fracture 0.00000 0.9306071 0.069393
## Death    0.00000 0.0000000 1.000000

1.3.5 Check the model goodness-of-fit

Model validation
prevalence.msm(multi.m2, times = seq(0,20,1), covariates = "mean", ci = "normal")
## $Observed
##    State 1 State 2 State 3 Total
## 0     5373       0       0  5373
## 1     5268      49      53  5370
## 2     5118     101     148  5367
## 3     4921     159     262  5342
## 4     4731     198     394  5323
## 5     4475     257     559  5291
## 6     4250     290     726  5266
## 7     4020     323     905  5248
## 8     3753     371    1102  5226
## 9     3495     392    1322  5209
## 10    3242     413    1533  5188
## 11    2984     419    1764  5167
## 12    2748     418    1978  5144
## 13    2480     414    2226  5120
## 14    2220     413    2455  5088
## 15    1961     399    2673  5033
## 16    1744     365    2871  4980
## 17    1528     346    3066  4940
## 18    1337     310    3230  4877
## 19    1181     265    3291  4737
## 20    1018     212    3292  4522
## 
## $Expected
## $Expected$estimates
##      Well Fracture    Death Total
## 0  5373.0    0.000    0.000  5373
## 1  5226.0   57.334   86.653  5370
## 2  5083.0  110.516  173.438  5367
## 3  4923.7  159.121  259.168  5342
## 4  4774.6  203.879  344.474  5323
## 5  4618.7  244.302  428.008  5291
## 6  4260.4  307.865  697.730  5266
## 7  3935.1  362.528  950.398  5248
## 8  3631.8  408.587 1185.651  5226
## 9  3355.0  447.570 1406.440  5209
## 10 3096.9  479.528 1611.581  5188
## 11 2858.6  505.488 1802.912  5167
## 12 2637.6  525.892 1980.533  5144
## 13 2433.1  541.414 2145.470  5120
## 14 2240.9  551.825 2295.242  5088
## 15 2054.5  555.909 2422.630  5033
## 16 1884.0  556.829 2539.135  4980
## 17 1732.1  556.281 2651.608  4940
## 18 1584.9  550.616 2741.525  4877
## 19 1426.7  534.101 2776.207  4737
## 20 1262.3  507.421 2752.327  4522
## 
## $Expected$ci
## , , 2.5%
## 
##         [,1]    [,2]     [,3] [,4]
##  [1,] 5373.0   0.000    0.000 5373
##  [2,] 5215.1  51.702   79.083 5370
##  [3,] 5061.9  99.710  158.530 5367
##  [4,] 4893.0 143.204  237.206 5342
##  [5,] 4734.9 183.219  315.927 5323
##  [6,] 4570.7 219.115  393.160 5291
##  [7,] 4213.3 284.004  663.601 5266
##  [8,] 3886.0 337.887  913.131 5248
##  [9,] 3581.9 382.569 1141.878 5226
## [10,] 3303.1 419.262 1357.583 5209
## [11,] 3041.1 448.918 1559.167 5188
## [12,] 2799.4 472.290 1745.416 5167
## [13,] 2576.1 490.427 1919.332 5144
## [14,] 2369.1 503.275 2081.167 5120
## [15,] 2175.8 512.157 2227.955 5088
## [16,] 1989.8 512.658 2352.849 5033
## [17,] 1820.6 511.819 2467.717 4980
## [18,] 1669.3 511.200 2578.596 4940
## [19,] 1523.3 505.227 2667.492 4877
## [20,] 1367.6 488.389 2702.586 4737
## [21,] 1206.8 462.175 2681.197 4522
## 
## , , 97.5%
## 
##         [,1]    [,2]     [,3] [,4]
##  [1,] 5373.0   0.000    0.000 5373
##  [2,] 5235.6  64.339   94.413 5370
##  [3,] 5101.6 124.042  188.778 5367
##  [4,] 4950.7 178.851  281.361 5342
##  [5,] 4809.6 229.484  373.338 5323
##  [6,] 4661.0 275.400  463.926 5291
##  [7,] 4302.7 337.960  734.380 5266
##  [8,] 3979.7 392.016  989.372 5248
##  [9,] 3679.8 438.461 1227.584 5226
## [10,] 3406.9 478.566 1452.528 5209
## [11,] 3152.8 512.715 1662.620 5188
## [12,] 2918.6 541.474 1857.465 5167
## [13,] 2701.1 564.304 2039.443 5144
## [14,] 2498.9 581.998 2208.432 5120
## [15,] 2306.7 595.035 2361.060 5088
## [16,] 2121.4 601.044 2491.266 5033
## [17,] 1951.5 603.343 2610.034 4980
## [18,] 1798.5 604.609 2724.109 4940
## [19,] 1650.0 599.487 2814.206 4877
## [20,] 1489.6 582.714 2848.005 4737
## [21,] 1321.8 555.839 2821.030 4522
## 
## 
## 
## $`Observed percentages`
##    State 1 State 2  State 3
## 0  100.000 0.00000  0.00000
## 1   98.101 0.91248  0.98696
## 2   95.361 1.88187  2.75759
## 3   92.119 2.97641  4.90453
## 4   88.878 3.71971  7.40184
## 5   84.578 4.85730 10.56511
## 6   80.706 5.50703 13.78656
## 7   76.601 6.15473 17.24466
## 8   71.814 7.09912 21.08687
## 9   67.095 7.52544 25.37915
## 10  62.490 7.96068 29.54896
## 11  57.751 8.10915 34.13973
## 12  53.421 8.12597 38.45257
## 13  48.438 8.08594 43.47656
## 14  43.632 8.11714 48.25079
## 15  38.963 7.92768 53.10948
## 16  35.020 7.32932 57.65060
## 17  30.931 7.00405 62.06478
## 18  27.414 6.35637 66.22924
## 19  24.931 5.59426 69.47435
## 20  22.512 4.68819 72.79965
## 
## $`Expected percentages`
## $`Expected percentages`$estimates
##       Well Fracture   Death
## 0  100.000   0.0000  0.0000
## 1   97.319   1.0677  1.6137
## 2   94.709   2.0592  3.2316
## 3   92.170   2.9787  4.8515
## 4   89.698   3.8301  6.4714
## 5   87.293   4.6173  8.0894
## 6   80.904   5.8463 13.2497
## 7   74.982   6.9079 18.1097
## 8   69.494   7.8184 22.6875
## 9   64.408   8.5922 27.0002
## 10  59.693   9.2430 31.0636
## 11  55.324   9.7830 34.8928
## 12  51.275  10.2234 38.5018
## 13  47.522  10.5745 41.9037
## 14  44.044  10.8456 45.1109
## 15  40.820  11.0453 48.1349
## 16  37.832  11.1813 50.9866
## 17  35.063  11.2608 53.6763
## 18  32.497  11.2901 56.2134
## 19  30.118  11.2751 58.6069
## 20  27.914  11.2212 60.8653
## 
## $`Expected percentages`$ci
## , , 2.5%
## 
##          [,1]     [,2]    [,3]
##  [1,] 100.000  0.00000  0.0000
##  [2,]  97.116  0.96279  1.4727
##  [3,]  94.315  1.85784  2.9538
##  [4,]  91.594  2.68071  4.4404
##  [5,]  88.953  3.44203  5.9351
##  [6,]  86.387  4.14127  7.4307
##  [7,]  80.010  5.39317 12.6016
##  [8,]  74.048  6.43840 17.3996
##  [9,]  68.540  7.32048 21.8499
## [10,]  63.412  8.04880 26.0623
## [11,]  58.619  8.65301 30.0533
## [12,]  54.179  9.14050 33.7801
## [13,]  50.080  9.53396 37.3121
## [14,]  46.272  9.82959 40.6478
## [15,]  42.763 10.06598 43.7884
## [16,]  39.536 10.18594 46.7484
## [17,]  36.558 10.27749 49.5526
## [18,]  33.792 10.34817 52.1983
## [19,]  31.235 10.35938 54.6953
## [20,]  28.872 10.31010 57.0527
## [21,]  26.687 10.22060 59.2923
## 
## , , 97.5%
## 
##          [,1]    [,2]    [,3]
##  [1,] 100.000  0.0000  0.0000
##  [2,]  97.496  1.1981  1.7581
##  [3,]  95.055  2.3112  3.5174
##  [4,]  92.675  3.3480  5.2670
##  [5,]  90.355  4.3112  7.0137
##  [6,]  88.093  5.2051  8.7682
##  [7,]  81.706  6.4178 13.9457
##  [8,]  75.833  7.4698 18.8524
##  [9,]  70.414  8.3900 23.4899
## [10,]  65.403  9.1873 27.8850
## [11,]  60.771  9.8827 32.0474
## [12,]  56.485 10.4795 35.9486
## [13,]  52.509 10.9701 39.6470
## [14,]  48.807 11.3672 43.1334
## [15,]  45.337 11.6949 46.4045
## [16,]  42.150 11.9421 49.4986
## [17,]  39.187 12.1153 52.4103
## [18,]  36.407 12.2391 55.1439
## [19,]  33.832 12.2921 57.7036
## [20,]  31.447 12.3013 60.1226
## [21,]  29.230 12.2919 62.3846
plot.prevalence.msm(multi.m2, mintime = 0, maxtime = 20, legend.pos = c(10, 80), col.obs = "gray", 
                    col.exp = "black")

Diagnostic plots for survival
par(mfrow= c(1,2))
plot.survfit.msm(multi.m2, from = 1, main = "No fracture", range = c(0,20), 
                 ci = "normal", col = "black", col.ci = "black", lty = 2, lwd = 2, 
                 col.surv ="grey", lty.surv = 2, lwd.surv = 1, xlab = "Time (years)")
legend(0.4,0.35, legend = c("Expected", "Expected (95% CI)", "Observed", "Observed (95% CI)"), 
       col = c("black", "black", "grey", "gray"), lty = c(2, 2, 1, 2), lwd = c(2, 1, 2, 1))
plot.survfit.msm(multi.m2, from = 2, main = "Initial fracture", range = c(0,20), 
                 ci = "normal", col = "black", col.ci = "black", lty = 2, lwd = 2, 
                 col.surv ="grey", lty.surv = 2, lwd.surv = 1, xlab = "Time (years)")

2. WOMEN

2.1 Description of study sample

women = subset(bmd, gender == "F" & race == "1:WHITE")
baseline.w = subset(women, state == 1 & time2event == 0)

table1(~ ageBase + fnbmdBase + TscoreBase + as.factor(fall.no) + fall.yesno + as.factor(fx50) + race + weight + height + BMI + as.factor(smoke) + as.factor(drink.n) + as.factor(physical) + as.factor(cvd.n) + as.factor(hypertension) + as.factor(copd) + as.factor(diabetes.n) + as.factor(cancer) + as.factor(parkinson) + as.factor(rheumatoid) + as.factor(renal) + as.factor(depression) + as.factor(anyfx) + as.factor(death) | fx_death, data = baseline.w)
Event-free
(N=1888)
Fractured, Alive
(N=1210)
Fractured, Dead
(N=1486)
No Fractured, Dead
(N=2779)
Overall
(N=7363)
ageBase
Mean (SD) 71.7 (4.00) 72.0 (4.20) 75.4 (5.29) 74.7 (5.33) 73.6 (5.08)
Median [Min, Max] 71.0 [67.0, 88.0] 71.0 [67.0, 90.0] 75.0 [67.0, 91.0] 74.0 [67.0, 90.0] 73.0 [67.0, 91.0]
fnbmdBase
Mean (SD) 0.676 (0.107) 0.638 (0.104) 0.604 (0.101) 0.658 (0.114) 0.648 (0.111)
Median [Min, Max] 0.666 [0.387, 1.15] 0.630 [0.302, 1.11] 0.597 [0.297, 1.13] 0.648 [0.277, 1.19] 0.639 [0.277, 1.19]
TscoreBase
Mean (SD) -1.54 (0.893) -1.85 (0.869) -2.13 (0.843) -1.69 (0.947) -1.76 (0.924)
Median [Min, Max] -1.62 [-3.94, 2.39] -1.92 [-4.65, 2.12] -2.20 [-4.69, 2.22] -1.77 [-4.86, 2.73] -1.84 [-4.86, 2.73]
as.factor(fall.no)
0 1418 (75.1%) 854 (70.6%) 1034 (69.6%) 2038 (73.3%) 5344 (72.6%)
1 333 (17.6%) 257 (21.2%) 286 (19.2%) 473 (17.0%) 1349 (18.3%)
2 97 (5.1%) 71 (5.9%) 107 (7.2%) 176 (6.3%) 451 (6.1%)
3+ 40 (2.1%) 28 (2.3%) 59 (4.0%) 92 (3.3%) 219 (3.0%)
fall.yesno
No 1418 (75.1%) 854 (70.6%) 1034 (69.6%) 2038 (73.3%) 5344 (72.6%)
Yes 470 (24.9%) 356 (29.4%) 452 (30.4%) 741 (26.7%) 2019 (27.4%)
as.factor(fx50)
0 1297 (68.7%) 654 (54.0%) 708 (47.6%) 1700 (61.2%) 4359 (59.2%)
1 591 (31.3%) 556 (46.0%) 778 (52.4%) 1079 (38.8%) 3004 (40.8%)
race
1:WHITE 1888 (100%) 1210 (100%) 1486 (100%) 2779 (100%) 7363 (100%)
weight
Mean (SD) 67.3 (11.8) 66.3 (11.7) 64.7 (11.6) 66.7 (12.8) 66.4 (12.1)
Median [Min, Max] 65.8 [41.7, 112] 65.0 [40.8, 112] 63.3 [40.2, 110] 65.3 [40.0, 112] 65.0 [40.0, 112]
Missing 11 (0.6%) 11 (0.9%) 20 (1.3%) 24 (0.9%) 66 (0.9%)
height
Mean (SD) 160 (5.66) 160 (5.65) 159 (6.01) 159 (5.94) 159 (5.85)
Median [Min, Max] 160 [142, 175] 160 [143, 176] 159 [142, 175] 159 [141, 175] 159 [141, 176]
Missing 8 (0.4%) 7 (0.6%) 16 (1.1%) 13 (0.5%) 44 (0.6%)
BMI
Mean (SD) 26.4 (4.36) 26.0 (4.37) 25.7 (4.35) 26.5 (4.69) 26.2 (4.50)
Median [Min, Max] 25.8 [16.0, 44.1] 25.4 [16.2, 48.8] 25.3 [16.9, 48.3] 25.9 [15.3, 44.7] 25.7 [15.3, 48.8]
Missing 11 (0.6%) 11 (0.9%) 20 (1.3%) 25 (0.9%) 67 (0.9%)
as.factor(smoke)
0 1182 (62.6%) 771 (63.7%) 915 (61.6%) 1614 (58.1%) 4482 (60.9%)
1 706 (37.4%) 439 (36.3%) 571 (38.4%) 1165 (41.9%) 2881 (39.1%)
as.factor(drink.n)
No 723 (38.3%) 480 (39.7%) 724 (48.7%) 1336 (48.1%) 3263 (44.3%)
Yes 1165 (61.7%) 730 (60.3%) 762 (51.3%) 1443 (51.9%) 4100 (55.7%)
as.factor(physical)
0 567 (30.0%) 355 (29.3%) 734 (49.4%) 1400 (50.4%) 3056 (41.5%)
1 1321 (70.0%) 855 (70.7%) 752 (50.6%) 1379 (49.6%) 4307 (58.5%)
as.factor(cvd.n)
No 1554 (82.3%) 1002 (82.8%) 1079 (72.6%) 1956 (70.4%) 5591 (75.9%)
Yes 334 (17.7%) 208 (17.2%) 407 (27.4%) 823 (29.6%) 1772 (24.1%)
as.factor(hypertension)
0 1326 (70.2%) 864 (71.4%) 846 (56.9%) 1555 (56.0%) 4591 (62.4%)
1 562 (29.8%) 346 (28.6%) 640 (43.1%) 1224 (44.0%) 2772 (37.6%)
as.factor(copd)
0 1628 (86.2%) 1007 (83.2%) 1190 (80.1%) 2289 (82.4%) 6114 (83.0%)
1 260 (13.8%) 203 (16.8%) 296 (19.9%) 490 (17.6%) 1249 (17.0%)
as.factor(diabetes.n)
No 1821 (96.5%) 1153 (95.3%) 1357 (91.3%) 2543 (91.5%) 6874 (93.4%)
Yes 67 (3.5%) 57 (4.7%) 129 (8.7%) 236 (8.5%) 489 (6.6%)
as.factor(cancer)
0 1546 (81.9%) 977 (80.7%) 1181 (79.5%) 2232 (80.3%) 5936 (80.6%)
1 342 (18.1%) 233 (19.3%) 305 (20.5%) 547 (19.7%) 1427 (19.4%)
as.factor(parkinson)
0 1882 (99.7%) 1204 (99.5%) 1470 (98.9%) 2764 (99.5%) 7320 (99.4%)
1 6 (0.3%) 6 (0.5%) 16 (1.1%) 15 (0.5%) 43 (0.6%)
as.factor(rheumatoid)
0 1804 (95.6%) 1146 (94.7%) 1393 (93.7%) 2629 (94.6%) 6972 (94.7%)
1 84 (4.4%) 64 (5.3%) 93 (6.3%) 150 (5.4%) 391 (5.3%)
as.factor(renal)
0 1867 (98.9%) 1203 (99.4%) 1463 (98.5%) 2741 (98.6%) 7274 (98.8%)
1 21 (1.1%) 7 (0.6%) 23 (1.5%) 38 (1.4%) 89 (1.2%)
as.factor(depression)
0 1765 (93.5%) 1117 (92.3%) 1352 (91.0%) 2562 (92.2%) 6796 (92.3%)
1 123 (6.5%) 93 (7.7%) 134 (9.0%) 217 (7.8%) 567 (7.7%)
as.factor(anyfx)
0 1888 (100%) 0 (0%) 0 (0%) 2779 (100%) 4667 (63.4%)
1 0 (0%) 1210 (100%) 1486 (100%) 0 (0%) 2696 (36.6%)
as.factor(death)
0 1888 (100%) 1210 (100%) 0 (0%) 0 (0%) 3098 (42.1%)
1 0 (0%) 0 (0%) 1486 (100%) 2779 (100%) 4265 (57.9%)
#createTable(compareGroups(fx_death ~ ageBase + fnbmdBase + TscoreBase + fall.no + fall.yesno + fx50 + race + weight + height + BMI + smoke + drink.n + physical + cvd.n + hypertension + copd + diabetes.n + cancer + parkinson + rheumatoid + renal + depression + anyfx + death, data = baseline.w))

2.2 Conventional Cox’s PH model

2.2.1 Age-adjusted model:
cox.w1 = coxph(Surv(time2end, death) ~ ageBase + TscoreBase.2, data = baseline.w)
summary(cox.w1)
## Call:
## coxph(formula = Surv(time2end, death) ~ ageBase + TscoreBase.2, 
##     data = baseline.w)
## 
##   n= 7363, number of events= 4265 
## 
##                 coef exp(coef) se(coef)     z Pr(>|z|)    
## ageBase      0.11955   1.12699  0.00308 38.83  < 2e-16 ***
## TscoreBase.2 0.07487   1.07774  0.01788  4.19  2.8e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##              exp(coef) exp(-coef) lower .95 upper .95
## ageBase           1.13      0.887      1.12      1.13
## TscoreBase.2      1.08      0.928      1.04      1.12
## 
## Concordance= 0.664  (se = 0.004 )
## Likelihood ratio test= 1545  on 2 df,   p=<2e-16
## Wald test            = 1718  on 2 df,   p=<2e-16
## Score (logrank) test = 1810  on 2 df,   p=<2e-16
2.2.2 Essential covariates
cox.w2 = coxph(Surv(time2end, death) ~ ageBase + TscoreBase.2 + fall.yesno + fx50 + cvd.n + hypertension + copd + diabetes.n + cancer, data = baseline.w)
summary(cox.w2)
## Call:
## coxph(formula = Surv(time2end, death) ~ ageBase + TscoreBase.2 + 
##     fall.yesno + fx50 + cvd.n + hypertension + copd + diabetes.n + 
##     cancer, data = baseline.w)
## 
##   n= 7363, number of events= 4265 
## 
##                   coef exp(coef) se(coef)     z Pr(>|z|)    
## ageBase        0.11275   1.11935  0.00316 35.72  < 2e-16 ***
## TscoreBase.2   0.09486   1.09951  0.01830  5.18  2.2e-07 ***
## fall.yesnoYes -0.02686   0.97349  0.03426 -0.78   0.4329    
## fx50           0.09454   1.09916  0.03184  2.97   0.0030 ** 
## cvd.nYes       0.32071   1.37811  0.03433  9.34  < 2e-16 ***
## hypertension   0.31682   1.37275  0.03137 10.10  < 2e-16 ***
## copd           0.12413   1.13216  0.03965  3.13   0.0017 ** 
## diabetes.nYes  0.58169   1.78905  0.05539 10.50  < 2e-16 ***
## cancer         0.02435   1.02465  0.03837  0.63   0.5257    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## ageBase           1.119      0.893      1.11      1.13
## TscoreBase.2      1.100      0.909      1.06      1.14
## fall.yesnoYes     0.973      1.027      0.91      1.04
## fx50              1.099      0.910      1.03      1.17
## cvd.nYes          1.378      0.726      1.29      1.47
## hypertension      1.373      0.728      1.29      1.46
## copd              1.132      0.883      1.05      1.22
## diabetes.nYes     1.789      0.559      1.60      1.99
## cancer            1.025      0.976      0.95      1.10
## 
## Concordance= 0.687  (se = 0.004 )
## Likelihood ratio test= 1877  on 9 df,   p=<2e-16
## Wald test            = 2065  on 9 df,   p=<2e-16
## Score (logrank) test = 2168  on 9 df,   p=<2e-16
2.2.3 Multivariable-adjusted model
cox.w3 = coxph(Surv(time2end, death) ~ ageBase + TscoreBase.2 + fall.yesno + fx50 + BMI + smoke + drink.n + 
                 physical + cvd.n + hypertension + copd + diabetes.n + cancer + renal + parkinson + 
                 depression, data = baseline.w)
summary(cox.w3)
## Call:
## coxph(formula = Surv(time2end, death) ~ ageBase + TscoreBase.2 + 
##     fall.yesno + fx50 + BMI + smoke + drink.n + physical + cvd.n + 
##     hypertension + copd + diabetes.n + cancer + renal + parkinson + 
##     depression, data = baseline.w)
## 
##   n= 7296, number of events= 4220 
##    (67 observations deleted due to missingness)
## 
##                   coef exp(coef) se(coef)      z Pr(>|z|)    
## ageBase        0.10628   1.11214  0.00327  32.51  < 2e-16 ***
## TscoreBase.2   0.06549   1.06768  0.01954   3.35  0.00080 ***
## fall.yesnoYes -0.00126   0.99874  0.03448  -0.04  0.97076    
## fx50           0.08041   1.08373  0.03218   2.50  0.01247 *  
## BMI           -0.01361   0.98648  0.00396  -3.44  0.00059 ***
## smoke          0.30303   1.35395  0.03264   9.28  < 2e-16 ***
## drink.nYes    -0.16704   0.84617  0.03238  -5.16  2.5e-07 ***
## physical      -0.59912   0.54930  0.03288 -18.22  < 2e-16 ***
## cvd.nYes       0.28535   1.33022  0.03462   8.24  < 2e-16 ***
## hypertension   0.29332   1.34088  0.03173   9.24  < 2e-16 ***
## copd           0.02906   1.02948  0.04045   0.72  0.47259    
## diabetes.nYes  0.47070   1.60112  0.05673   8.30  < 2e-16 ***
## cancer         0.06810   1.07047  0.03869   1.76  0.07842 .  
## renal          0.36107   1.43486  0.13048   2.77  0.00565 ** 
## parkinson      0.20283   1.22487  0.18680   1.09  0.27755    
## depression     0.08250   1.08600  0.05664   1.46  0.14519    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## ageBase           1.112      0.899     1.105     1.119
## TscoreBase.2      1.068      0.937     1.028     1.109
## fall.yesnoYes     0.999      1.001     0.933     1.069
## fx50              1.084      0.923     1.017     1.154
## BMI               0.986      1.014     0.979     0.994
## smoke             1.354      0.739     1.270     1.443
## drink.nYes        0.846      1.182     0.794     0.902
## physical          0.549      1.821     0.515     0.586
## cvd.nYes          1.330      0.752     1.243     1.424
## hypertension      1.341      0.746     1.260     1.427
## copd              1.029      0.971     0.951     1.114
## diabetes.nYes     1.601      0.625     1.433     1.789
## cancer            1.070      0.934     0.992     1.155
## renal             1.435      0.697     1.111     1.853
## parkinson         1.225      0.816     0.849     1.766
## depression        1.086      0.921     0.972     1.213
## 
## Concordance= 0.716  (se = 0.004 )
## Likelihood ratio test= 2329  on 16 df,   p=<2e-16
## Wald test            = 2502  on 16 df,   p=<2e-16
## Score (logrank) test = 2653  on 16 df,   p=<2e-16

2.3 Multistate regression

2.3.1 Data preparation

qmatrix.n = crudeinits.msm(state ~ time2event, subject = ID, data = women, qmatrix = Q.m2)
qmatrix.n
##               Well  Fracture    Death
## Well     -0.069616  0.034280 0.035335
## Fracture  0.000000 -0.081983 0.081983
## Death     0.000000  0.000000 0.000000
statetable.msm(state, ID, data = women)
##     to
## from    1    2    3
##    1 1888 2696 2779
##    2    0 1210 1486

2.3.2 Age-adjusted model:

Age and FNBMD at baseline:
age.w1= msm(state ~ time2event, subject = ID, data = women, qmatrix = qmatrix.n, gen.inits = TRUE, exacttimes = TRUE, 
           method = "BFGS", control = list(fnscale = 4000, maxit = 10000), pci = 5, covariates =~ ageBase + TscoreBase.2)  
age.w1
## 
## Call:
## msm(formula = state ~ time2event, subject = ID, data = women,     qmatrix = qmatrix.n, gen.inits = TRUE, covariates = ~ageBase +         TscoreBase.2, exacttimes = TRUE, pci = 5, method = "BFGS",     control = list(fnscale = 4000, maxit = 10000))
## 
## Maximum likelihood estimates
## Baselines are with covariates set to their means
## 
## Transition intensities with hazard ratios for each covariate
##                     Baseline                     ageBase            
## Well - Well         -0.06128 (-0.06320,-0.05942)                    
## Well - Fracture      0.03321 ( 0.03189, 0.03458) 1.040 (1.032,1.048)
## Well - Death         0.02807 ( 0.02677, 0.02944) 1.105 (1.097,1.113)
## Fracture - Fracture -0.05900 (-0.06485,-0.05369)                    
## Fracture - Death     0.05900 ( 0.05369, 0.06485) 1.096 (1.085,1.107)
##                     TscoreBase.2          timeperiod[5,Inf)  
## Well - Well                                                  
## Well - Fracture     1.5795 (1.5067,1.656) 1.308 (1.209,1.416)
## Well - Death        0.9858 (0.9452,1.028) 3.093 (2.821,3.390)
## Fracture - Fracture                                          
## Fracture - Death    1.0820 (1.0134,1.155) 1.899 (1.592,2.264)
## 
## -2 * log-likelihood:  55919
hazard.msm(age.w1)
## $ageBase
##                      HR      L      U
## Well - Fracture  1.0402 1.0323 1.0483
## Well - Death     1.1051 1.0972 1.1131
## Fracture - Death 1.0961 1.0855 1.1069
## 
## $TscoreBase.2
##                       HR       L      U
## Well - Fracture  1.57949 1.50667 1.6558
## Well - Death     0.98579 0.94516 1.0282
## Fracture - Death 1.08203 1.01342 1.1553
## 
## $`timeperiod[5,Inf)`
##                      HR      L      U
## Well - Fracture  1.3084 1.2091 1.4159
## Well - Death     3.0926 2.8209 3.3904
## Fracture - Death 1.8986 1.5920 2.2641
Age and FNBMD at baseline and at fracture time:
age.w2= msm(state ~ time2event, subject = ID, data = women, qmatrix = qmatrix.n, gen.inits = TRUE, exacttimes = TRUE, 
           method = "BFGS", control = list(fnscale = 4000, maxit = 10000), pci = 5, covariates =~ age + Tscore.2)  
age.w2
## 
## Call:
## msm(formula = state ~ time2event, subject = ID, data = women,     qmatrix = qmatrix.n, gen.inits = TRUE, covariates = ~age +         Tscore.2, exacttimes = TRUE, pci = 5, method = "BFGS",     control = list(fnscale = 4000, maxit = 10000))
## 
## Maximum likelihood estimates
## Baselines are with covariates set to their means
## 
## Transition intensities with hazard ratios for each covariate
##                     Baseline                     age                
## Well - Well         -0.06683 (-0.06892,-0.06480)                    
## Well - Fracture      0.03494 ( 0.03353, 0.03642) 1.041 (1.033,1.049)
## Well - Death         0.03189 ( 0.03045, 0.03339) 1.106 (1.098,1.114)
## Fracture - Fracture -0.04280 (-0.04751,-0.03855)                    
## Fracture - Death     0.04280 ( 0.03855, 0.04751) 1.105 (1.096,1.114)
##                     Tscore.2              timeperiod[5,Inf)  
## Well - Well                                                  
## Well - Fracture     1.5757 (1.5030,1.652) 1.312 (1.212,1.420)
## Well - Death        0.9811 (0.9406,1.023) 3.132 (2.855,3.435)
## Fracture - Fracture                                          
## Fracture - Death    1.1221 (1.0527,1.196) 1.196 (1.000,1.429)
## 
## -2 * log-likelihood:  55414
hazard.msm(age.w2)
## $age
##                      HR      L      U
## Well - Fracture  1.0405 1.0325 1.0486
## Well - Death     1.1061 1.0981 1.1141
## Fracture - Death 1.1047 1.0956 1.1139
## 
## $Tscore.2
##                       HR       L      U
## Well - Fracture  1.57570 1.50295 1.6520
## Well - Death     0.98111 0.94061 1.0234
## Fracture - Death 1.12210 1.05275 1.1960
## 
## $`timeperiod[5,Inf)`
##                      HR      L      U
## Well - Fracture  1.3118 1.2119 1.4199
## Well - Death     3.1317 2.8548 3.4354
## Fracture - Death 1.1955 1.0003 1.4288

2.3.3 Essential covariates:

Age and FNBMD at baseline:
multi.we1= msm(state ~ time2event, subject = ID, data = women, qmatrix = qmatrix.n, gen.inits = TRUE, exacttimes = TRUE, 
           method = "BFGS", control = list(fnscale = 4000, maxit = 10000), pci = 5, covariates =~ ageBase + TscoreBase.2 +
             fall.yesno + fx50 + cvd.n + hypertension + copd + diabetes.n + cancer)  
multi.we1
## 
## Call:
## msm(formula = state ~ time2event, subject = ID, data = women,     qmatrix = qmatrix.n, gen.inits = TRUE, covariates = ~ageBase +         TscoreBase.2 + fall.yesno + fx50 + cvd.n + hypertension +         copd + diabetes.n + cancer, exacttimes = TRUE, pci = 5,     method = "BFGS", control = list(fnscale = 4000, maxit = 10000))
## 
## Maximum likelihood estimates
## Baselines are with covariates set to their means
## 
## Transition intensities with hazard ratios for each covariate
##                     Baseline                     ageBase            
## Well - Well         -0.06070 (-0.06262,-0.05883)                    
## Well - Fracture      0.03306 ( 0.03174, 0.03444) 1.036 (1.028,1.044)
## Well - Death         0.02763 ( 0.02634, 0.02900) 1.097 (1.089,1.105)
## Fracture - Fracture -0.05819 (-0.06402,-0.05289)                    
## Fracture - Death     0.05819 ( 0.05289, 0.06402) 1.090 (1.079,1.101)
##                     TscoreBase.2         fall.yesnoYes        
## Well - Well                                                   
## Well - Fracture     1.540 (1.4678,1.616) 1.1510 (1.0596,1.250)
## Well - Death        1.009 (0.9661,1.053) 0.9697 (0.8912,1.055)
## Fracture - Fracture                                           
## Fracture - Death    1.107 (1.0361,1.183) 0.9416 (0.8425,1.052)
##                     fx50                  cvd.nYes             
## Well - Well                                                    
## Well - Fracture     1.4806 (1.3703,1.600) 0.9945 (0.9079,1.089)
## Well - Death        1.0771 (0.9960,1.165) 1.3773 (1.2682,1.496)
## Fracture - Fracture                                            
## Fracture - Death    0.9459 (0.8518,1.050) 1.2497 (1.1121,1.404)
##                     hypertension         copd                
## Well - Well                                                  
## Well - Fracture     1.072 (0.9896,1.160) 1.133 (1.0277,1.249)
## Well - Death        1.330 (1.2329,1.435) 1.132 (1.0265,1.248)
## Fracture - Fracture                                          
## Fracture - Death    1.310 (1.1787,1.456) 1.096 (0.9638,1.247)
##                     diabetes.nYes       cancer               
## Well - Well                                                  
## Well - Fracture     1.473 (1.267,1.712) 1.0472 (0.9527,1.151)
## Well - Death        1.601 (1.399,1.832) 1.0076 (0.9175,1.107)
## Fracture - Fracture                                          
## Fracture - Death    1.702 (1.417,2.044) 0.9854 (0.8685,1.118)
##                     timeperiod[5,Inf)  
## Well - Well                            
## Well - Fracture     1.343 (1.241,1.454)
## Well - Death        3.197 (2.916,3.506)
## Fracture - Fracture                    
## Fracture - Death    1.952 (1.636,2.328)
## 
## -2 * log-likelihood:  55507
hazard.msm(multi.we1)
## $ageBase
##                      HR      L      U
## Well - Fracture  1.0362 1.0280 1.0444
## Well - Death     1.0973 1.0892 1.1055
## Fracture - Death 1.0900 1.0790 1.1012
## 
## $TscoreBase.2
##                      HR      L      U
## Well - Fracture  1.5401 1.4678 1.6159
## Well - Death     1.0086 0.9661 1.0530
## Fracture - Death 1.1072 1.0361 1.1832
## 
## $fall.yesnoYes
##                       HR       L      U
## Well - Fracture  1.15105 1.05957 1.2504
## Well - Death     0.96975 0.89121 1.0552
## Fracture - Death 0.94160 0.84247 1.0524
## 
## $fx50
##                       HR       L      U
## Well - Fracture  1.48061 1.37026 1.5998
## Well - Death     1.07711 0.99604 1.1648
## Fracture - Death 0.94592 0.85177 1.0505
## 
## $cvd.nYes
##                       HR       L      U
## Well - Fracture  0.99449 0.90795 1.0893
## Well - Death     1.37731 1.26818 1.4958
## Fracture - Death 1.24974 1.11207 1.4045
## 
## $hypertension
##                      HR       L      U
## Well - Fracture  1.0715 0.98961 1.1602
## Well - Death     1.3301 1.23293 1.4350
## Fracture - Death 1.3101 1.17870 1.4562
## 
## $copd
##                      HR       L      U
## Well - Fracture  1.1331 1.02772 1.2492
## Well - Death     1.1321 1.02654 1.2485
## Fracture - Death 1.0963 0.96376 1.2470
## 
## $diabetes.nYes
##                      HR      L      U
## Well - Fracture  1.4730 1.2674 1.7119
## Well - Death     1.6013 1.3993 1.8324
## Fracture - Death 1.7019 1.4173 2.0438
## 
## $cancer
##                       HR       L      U
## Well - Fracture  1.04719 0.95266 1.1511
## Well - Death     1.00762 0.91750 1.1066
## Fracture - Death 0.98537 0.86852 1.1179
## 
## $`timeperiod[5,Inf)`
##                      HR      L      U
## Well - Fracture  1.3435 1.2412 1.4542
## Well - Death     3.1975 2.9159 3.5062
## Fracture - Death 1.9517 1.6360 2.3284
Age and FNBMD at fracture time:
multi.we2= msm(state ~ time2event, subject = ID, data = women, qmatrix = qmatrix.n, gen.inits = TRUE, exacttimes = TRUE, 
           method = "BFGS", control = list(fnscale = 4000, maxit = 10000), pci = 5, covariates =~ age + Tscore.2 +
             fall.yesno + fx50 + cvd.n + hypertension + copd + diabetes.n + cancer)  
multi.we2
## 
## Call:
## msm(formula = state ~ time2event, subject = ID, data = women,     qmatrix = qmatrix.n, gen.inits = TRUE, covariates = ~age +         Tscore.2 + fall.yesno + fx50 + cvd.n + hypertension +         copd + diabetes.n + cancer, exacttimes = TRUE, pci = 5,     method = "BFGS", control = list(fnscale = 4000, maxit = 10000))
## 
## Maximum likelihood estimates
## Baselines are with covariates set to their means
## 
## Transition intensities with hazard ratios for each covariate
##                     Baseline                     age                
## Well - Well         -0.06570 (-0.06779,-0.06368)                    
## Well - Fracture      0.03461 ( 0.03319, 0.03609) 1.036 (1.028,1.045)
## Well - Death         0.03109 ( 0.02967, 0.03259) 1.098 (1.090,1.107)
## Fracture - Fracture -0.04158 (-0.04624,-0.03739)                    
## Fracture - Death     0.04158 ( 0.03739, 0.04624) 1.103 (1.093,1.112)
##                     Tscore.2             fall.yesnoYes        
## Well - Well                                                   
## Well - Fracture     1.537 (1.4643,1.612) 1.1442 (1.0529,1.243)
## Well - Death        1.004 (0.9616,1.048) 0.9707 (0.8919,1.057)
## Fracture - Fracture                                           
## Fracture - Death    1.129 (1.0582,1.204) 0.9511 (0.8504,1.064)
##                     fx50                 cvd.nYes             
## Well - Well                                                   
## Well - Fracture     1.481 (1.3708,1.601) 0.9931 (0.9064,1.088)
## Well - Death        1.074 (0.9931,1.162) 1.3763 (1.2669,1.495)
## Fracture - Fracture                                           
## Fracture - Death    1.038 (0.9350,1.153) 1.2836 (1.1430,1.442)
##                     hypertension         copd                
## Well - Well                                                  
## Well - Fracture     1.072 (0.9895,1.161) 1.131 (1.0260,1.248)
## Well - Death        1.320 (1.2230,1.424) 1.140 (1.0335,1.257)
## Fracture - Fracture                                          
## Fracture - Death    1.331 (1.1986,1.479) 1.103 (0.9696,1.255)
##                     diabetes.nYes       cancer              
## Well - Well                                                 
## Well - Fracture     1.467 (1.262,1.706) 1.048 (0.9536,1.153)
## Well - Death        1.611 (1.407,1.843) 1.014 (0.9227,1.113)
## Fracture - Fracture                                         
## Fracture - Death    1.773 (1.476,2.131) 1.052 (0.9266,1.194)
##                     timeperiod[5,Inf)  
## Well - Well                            
## Well - Fracture     1.347 (1.244,1.458)
## Well - Death        3.237 (2.950,3.552)
## Fracture - Fracture                    
## Fracture - Death    1.277 (1.067,1.527)
## 
## -2 * log-likelihood:  54994
hazard.msm(multi.we2)
## $age
##                      HR      L      U
## Well - Fracture  1.0364 1.0283 1.0447
## Well - Death     1.0985 1.0903 1.1066
## Fracture - Death 1.1025 1.0933 1.1119
## 
## $Tscore.2
##                      HR       L      U
## Well - Fracture  1.5366 1.46433 1.6124
## Well - Death     1.0040 0.96164 1.0483
## Fracture - Death 1.1290 1.05821 1.2044
## 
## $fall.yesnoYes
##                       HR       L      U
## Well - Fracture  1.14419 1.05295 1.2433
## Well - Death     0.97071 0.89188 1.0565
## Fracture - Death 0.95107 0.85039 1.0637
## 
## $fx50
##                      HR       L      U
## Well - Fracture  1.4814 1.37080 1.6010
## Well - Death     1.0743 0.99313 1.1620
## Fracture - Death 1.0385 0.93496 1.1534
## 
## $cvd.nYes
##                       HR       L      U
## Well - Fracture  0.99311 0.90644 1.0881
## Well - Death     1.37632 1.26693 1.4952
## Fracture - Death 1.28364 1.14296 1.4416
## 
## $hypertension
##                      HR       L      U
## Well - Fracture  1.0716 0.98952 1.1605
## Well - Death     1.3198 1.22301 1.4243
## Fracture - Death 1.3315 1.19859 1.4791
## 
## $copd
##                      HR       L      U
## Well - Fracture  1.1314 1.02598 1.2476
## Well - Death     1.1399 1.03355 1.2572
## Fracture - Death 1.1030 0.96959 1.2547
## 
## $diabetes.nYes
##                      HR      L      U
## Well - Fracture  1.4672 1.2619 1.7060
## Well - Death     1.6106 1.4074 1.8431
## Fracture - Death 1.7731 1.4756 2.1307
## 
## $cancer
##                      HR       L      U
## Well - Fracture  1.0484 0.95360 1.1526
## Well - Death     1.0135 0.92274 1.1132
## Fracture - Death 1.0519 0.92656 1.1943
## 
## $`timeperiod[5,Inf)`
##                      HR      L      U
## Well - Fracture  1.3467 1.2439 1.4580
## Well - Death     3.2372 2.9502 3.5520
## Fracture - Death 1.2766 1.0671 1.5272

2.3.4 Multivariable-adjusted model:

All covariates at baseline:
multi.w1= msm(state ~ time2event, subject = ID, data = women, qmatrix = qmatrix.n, gen.inits = TRUE, exacttimes = TRUE, 
           method = "BFGS", control = list(fnscale = 4000, maxit = 10000), pci = 5, covariates =~ ageBase + TscoreBase.2 +
             fall.yesno + fx50 + BMI + smoke + drink.n + physical + cvd.n + hypertension + copd + diabetes.n + cancer + 
             renal + parkinson + depression)  
multi.w1
## 
## Call:
## msm(formula = state ~ time2event, subject = ID, data = women,     qmatrix = qmatrix.n, gen.inits = TRUE, covariates = ~ageBase +         TscoreBase.2 + fall.yesno + fx50 + BMI + smoke + drink.n +         physical + cvd.n + hypertension + copd + diabetes.n +         cancer + renal + parkinson + depression, exacttimes = TRUE,     pci = 5, method = "BFGS", control = list(fnscale = 4000,         maxit = 10000))
## 
## Maximum likelihood estimates
## Baselines are with covariates set to their means
## 
## Transition intensities with hazard ratios for each covariate
##                     Baseline                     ageBase            
## Well - Well         -0.05973 (-0.06165,-0.05787)                    
## Well - Fracture      0.03295 ( 0.03163, 0.03433) 1.032 (1.024,1.041)
## Well - Death         0.02678 ( 0.02549, 0.02813) 1.088 (1.080,1.096)
## Fracture - Fracture -0.05636 (-0.06212,-0.05113)                    
## Fracture - Death     0.05636 ( 0.05113, 0.06212) 1.087 (1.076,1.099)
##                     TscoreBase.2          fall.yesnoYes       
## Well - Well                                                   
## Well - Fracture     1.5290 (1.4519,1.610) 1.146 (1.0537,1.245)
## Well - Death        0.9881 (0.9437,1.035) 0.999 (0.9177,1.087)
## Fracture - Fracture                                           
## Fracture - Death    1.0547 (0.9815,1.133) 0.953 (0.8517,1.066)
##                     fx50                  BMI                   
## Well - Well                                                     
## Well - Fracture     1.4770 (1.3658,1.597) 0.9957 (0.9861,1.0055)
## Well - Death        1.0615 (0.9809,1.149) 0.9876 (0.9784,0.9969)
## Fracture - Fracture                                             
## Fracture - Death    0.9282 (0.8344,1.032) 0.9898 (0.9765,1.0034)
##                     smoke                 drink.nYes            
## Well - Well                                                     
## Well - Fracture     0.9911 (0.9136,1.075) 0.9731 (0.8984,1.0540)
## Well - Death        1.3373 (1.2361,1.447) 0.8724 (0.8067,0.9434)
## Fracture - Fracture                                             
## Fracture - Death    1.2351 (1.1056,1.380) 0.8260 (0.7404,0.9216)
##                     physical               cvd.nYes             
## Well - Well                                                     
## Well - Fracture     0.8309 (0.7654,0.9019) 0.9826 (0.8964,1.077)
## Well - Death        0.5633 (0.5203,0.6098) 1.3234 (1.2178,1.438)
## Fracture - Fracture                                             
## Fracture - Death    0.6508 (0.5834,0.7260) 1.2125 (1.0773,1.365)
##                     hypertension        copd                
## Well - Well                                                 
## Well - Fracture     1.057 (0.975,1.146) 1.108 (1.0030,1.224)
## Well - Death        1.316 (1.219,1.421) 1.043 (0.9438,1.152)
## Fracture - Fracture                                         
## Fracture - Death    1.251 (1.124,1.392) 1.028 (0.9015,1.173)
##                     diabetes.nYes       cancer              
## Well - Well                                                 
## Well - Fracture     1.430 (1.227,1.667) 1.063 (0.9666,1.170)
## Well - Death        1.427 (1.243,1.639) 1.057 (0.9621,1.162)
## Fracture - Fracture                                         
## Fracture - Death    1.538 (1.274,1.856) 1.008 (0.8871,1.145)
##                     renal                parkinson            
## Well - Well                                                   
## Well - Fracture     1.048 (0.7304,1.504) 2.0886 (1.3579,3.213)
## Well - Death        1.222 (0.8820,1.693) 1.2961 (0.7630,2.202)
## Fracture - Fracture                                           
## Fracture - Death    1.556 (1.0275,2.357) 0.9516 (0.5705,1.587)
##                     depression           timeperiod[5,Inf)  
## Well - Well                                                 
## Well - Fracture     1.115 (0.9707,1.280) 1.371 (1.266,1.485)
## Well - Death        1.030 (0.8949,1.185) 3.368 (3.069,3.696)
## Fracture - Fracture                                         
## Fracture - Death    1.026 (0.8553,1.231) 2.109 (1.762,2.523)
## 
## -2 * log-likelihood:  54539
hazard.msm(multi.w1)
## $ageBase
##                      HR      L      U
## Well - Fracture  1.0322 1.0238 1.0406
## Well - Death     1.0880 1.0796 1.0964
## Fracture - Death 1.0874 1.0761 1.0988
## 
## $TscoreBase.2
##                       HR       L      U
## Well - Fracture  1.52904 1.45190 1.6103
## Well - Death     0.98813 0.94366 1.0347
## Fracture - Death 1.05467 0.98146 1.1333
## 
## $fall.yesnoYes
##                       HR       L      U
## Well - Fracture  1.14555 1.05372 1.2454
## Well - Death     0.99897 0.91770 1.0874
## Fracture - Death 0.95299 0.85169 1.0663
## 
## $fx50
##                       HR       L      U
## Well - Fracture  1.47704 1.36578 1.5974
## Well - Death     1.06153 0.98089 1.1488
## Fracture - Death 0.92818 0.83445 1.0324
## 
## $BMI
##                       HR       L       U
## Well - Fracture  0.99574 0.98608 1.00548
## Well - Death     0.98761 0.97838 0.99692
## Fracture - Death 0.98985 0.97645 1.00342
## 
## $smoke
##                       HR       L      U
## Well - Fracture  0.99114 0.91364 1.0752
## Well - Death     1.33728 1.23610 1.4467
## Fracture - Death 1.23508 1.10562 1.3797
## 
## $drink.nYes
##                       HR       L       U
## Well - Fracture  0.97312 0.89841 1.05404
## Well - Death     0.87238 0.80669 0.94342
## Fracture - Death 0.82604 0.74037 0.92163
## 
## $physical
##                       HR       L       U
## Well - Fracture  0.83087 0.76541 0.90193
## Well - Death     0.56327 0.52031 0.60977
## Fracture - Death 0.65081 0.58341 0.72600
## 
## $cvd.nYes
##                       HR       L      U
## Well - Fracture  0.98259 0.89639 1.0771
## Well - Death     1.32340 1.21780 1.4382
## Fracture - Death 1.21250 1.07733 1.3646
## 
## $hypertension
##                      HR       L      U
## Well - Fracture  1.0568 0.97503 1.1455
## Well - Death     1.3159 1.21869 1.4210
## Fracture - Death 1.2508 1.12379 1.3923
## 
## $copd
##                      HR       L      U
## Well - Fracture  1.1081 1.00300 1.2242
## Well - Death     1.0427 0.94384 1.1519
## Fracture - Death 1.0283 0.90147 1.1730
## 
## $diabetes.nYes
##                      HR      L      U
## Well - Fracture  1.4301 1.2267 1.6671
## Well - Death     1.4271 1.2429 1.6386
## Fracture - Death 1.5378 1.2739 1.8565
## 
## $cancer
##                      HR       L      U
## Well - Fracture  1.0635 0.96663 1.1700
## Well - Death     1.0574 0.96213 1.1621
## Fracture - Death 1.0078 0.88709 1.1450
## 
## $renal
##                      HR       L      U
## Well - Fracture  1.0480 0.73036 1.5037
## Well - Death     1.2220 0.88195 1.6931
## Fracture - Death 1.5562 1.02753 2.3569
## 
## $parkinson
##                       HR       L      U
## Well - Fracture  2.08862 1.35787 3.2126
## Well - Death     1.29614 0.76302 2.2018
## Fracture - Death 0.95156 0.57052 1.5871
## 
## $depression
##                      HR       L      U
## Well - Fracture  1.1149 0.97070 1.2805
## Well - Death     1.0299 0.89494 1.1853
## Fracture - Death 1.0260 0.85532 1.2308
## 
## $`timeperiod[5,Inf)`
##                      HR      L      U
## Well - Fracture  1.3710 1.2656 1.4852
## Well - Death     3.3683 3.0693 3.6963
## Fracture - Death 2.1086 1.7624 2.5227
Age and FNBMD at fracture time (other covariates at baseline):
multi.w2= msm(state ~ time2event, subject = ID, data = women, qmatrix = qmatrix.n, gen.inits = TRUE, exacttimes = TRUE, 
           method = "BFGS", control = list(fnscale = 4000, maxit = 10000), pci = 5, covariates =~ age + Tscore.2 + fall.yesno
           + fx50 + BMI + smoke + drink.n + physical + cvd.n + hypertension + copd + diabetes.n + cancer + renal + parkinson 
           + depression)  
multi.w2
## 
## Call:
## msm(formula = state ~ time2event, subject = ID, data = women,     qmatrix = qmatrix.n, gen.inits = TRUE, covariates = ~age +         Tscore.2 + fall.yesno + fx50 + BMI + smoke + drink.n +         physical + cvd.n + hypertension + copd + diabetes.n +         cancer + renal + parkinson + depression, exacttimes = TRUE,     pci = 5, method = "BFGS", control = list(fnscale = 4000,         maxit = 10000))
## 
## Maximum likelihood estimates
## Baselines are with covariates set to their means
## 
## Transition intensities with hazard ratios for each covariate
##                     Baseline                     age                
## Well - Well         -0.06416 (-0.06624,-0.06215)                    
## Well - Fracture      0.03431 ( 0.03289, 0.03579) 1.032 (1.024,1.041)
## Well - Death         0.02986 ( 0.02845, 0.03134) 1.089 (1.081,1.098)
## Fracture - Fracture -0.03959 (-0.04414,-0.03551)                    
## Fracture - Death     0.03959 ( 0.03551, 0.04414) 1.104 (1.094,1.113)
##                     Tscore.2              fall.yesnoYes        
## Well - Well                                                    
## Well - Fracture     1.5240 (1.4470,1.605) 1.1391 (1.0475,1.239)
## Well - Death        0.9864 (0.9419,1.033) 0.9992 (0.9177,1.088)
## Fracture - Fracture                                            
## Fracture - Death    1.0633 (0.9911,1.141) 0.9593 (0.8566,1.074)
##                     fx50                 BMI                   
## Well - Well                                                    
## Well - Fracture     1.478 (1.3669,1.599) 0.9954 (0.9857,1.0051)
## Well - Death        1.058 (0.9773,1.145) 0.9883 (0.9791,0.9976)
## Fracture - Fracture                                            
## Fracture - Death    1.021 (0.9178,1.136) 0.9862 (0.9728,0.9997)
##                     smoke                drink.nYes            
## Well - Well                                                    
## Well - Fracture     0.995 (0.9171,1.080) 0.9724 (0.8976,1.0534)
## Well - Death        1.347 (1.2447,1.457) 0.8737 (0.8077,0.9451)
## Fracture - Fracture                                            
## Fracture - Death    1.289 (1.1538,1.441) 0.8055 (0.7216,0.8992)
##                     physical               cvd.nYes            
## Well - Well                                                    
## Well - Fracture     0.8320 (0.7663,0.9033) 0.9813 (0.895,1.076)
## Well - Death        0.5685 (0.5251,0.6156) 1.3217 (1.216,1.437)
## Fracture - Fracture                                            
## Fracture - Death    0.6352 (0.5697,0.7081) 1.2404 (1.103,1.395)
##                     hypertension         copd                
## Well - Well                                                  
## Well - Fracture     1.058 (0.9757,1.147) 1.106 (1.0008,1.222)
## Well - Death        1.307 (1.2105,1.412) 1.048 (0.9485,1.158)
## Fracture - Fracture                                          
## Fracture - Death    1.268 (1.1396,1.411) 1.031 (0.9036,1.176)
##                     diabetes.nYes       cancer              
## Well - Well                                                 
## Well - Fracture     1.426 (1.223,1.663) 1.065 (0.9675,1.171)
## Well - Death        1.435 (1.250,1.648) 1.062 (0.9664,1.168)
## Fracture - Fracture                                         
## Fracture - Death    1.572 (1.300,1.900) 1.079 (0.9494,1.227)
##                     renal                parkinson           
## Well - Well                                                  
## Well - Fracture     1.051 (0.7323,1.508) 2.095 (1.3621,3.222)
## Well - Death        1.229 (0.8874,1.703) 1.306 (0.7696,2.216)
## Fracture - Fracture                                          
## Fracture - Death    1.763 (1.1628,2.673) 1.098 (0.6580,1.833)
##                     depression           timeperiod[5,Inf)  
## Well - Well                                                 
## Well - Fracture     1.112 (0.9682,1.278) 1.374 (1.268,1.489)
## Well - Death        1.035 (0.8989,1.191) 3.395 (3.092,3.728)
## Fracture - Fracture                                         
## Fracture - Death    1.072 (0.8925,1.287) 1.400 (1.167,1.681)
## 
## -2 * log-likelihood:  54034
hazard.msm(multi.w2)
## $age
##                      HR      L      U
## Well - Fracture  1.0324 1.0240 1.0409
## Well - Death     1.0892 1.0808 1.0977
## Fracture - Death 1.1037 1.0942 1.1132
## 
## $Tscore.2
##                       HR      L      U
## Well - Fracture  1.52397 1.4470 1.6051
## Well - Death     0.98638 0.9419 1.0330
## Fracture - Death 1.06326 0.9911 1.1407
## 
## $fall.yesnoYes
##                       HR       L      U
## Well - Fracture  1.13912 1.04751 1.2388
## Well - Death     0.99918 0.91768 1.0879
## Fracture - Death 0.95927 0.85661 1.0742
## 
## $fx50
##                      HR       L      U
## Well - Fracture  1.4785 1.36687 1.5992
## Well - Death     1.0579 0.97728 1.1451
## Fracture - Death 1.0212 0.91781 1.1362
## 
## $BMI
##                       HR       L       U
## Well - Fracture  0.99538 0.98571 1.00514
## Well - Death     0.98831 0.97906 0.99764
## Fracture - Death 0.98618 0.97283 0.99972
## 
## $smoke
##                       HR       L      U
## Well - Fracture  0.99497 0.91707 1.0795
## Well - Death     1.34685 1.24472 1.4574
## Fracture - Death 1.28944 1.15382 1.4410
## 
## $drink.nYes
##                       HR       L       U
## Well - Fracture  0.97239 0.89758 1.05343
## Well - Death     0.87372 0.80775 0.94507
## Fracture - Death 0.80553 0.72161 0.89922
## 
## $physical
##                       HR       L       U
## Well - Fracture  0.83195 0.76625 0.90329
## Well - Death     0.56852 0.52505 0.61558
## Fracture - Death 0.63516 0.56972 0.70811
## 
## $cvd.nYes
##                       HR       L      U
## Well - Fracture  0.98135 0.89501 1.0760
## Well - Death     1.32169 1.21591 1.4367
## Fracture - Death 1.24037 1.10258 1.3954
## 
## $hypertension
##                      HR       L      U
## Well - Fracture  1.0578 0.97571 1.1468
## Well - Death     1.3074 1.21046 1.4120
## Fracture - Death 1.2680 1.13959 1.4108
## 
## $copd
##                      HR       L      U
## Well - Fracture  1.1059 1.00079 1.2221
## Well - Death     1.0479 0.94847 1.1577
## Fracture - Death 1.0307 0.90359 1.1758
## 
## $diabetes.nYes
##                      HR      L      U
## Well - Fracture  1.4258 1.2225 1.6628
## Well - Death     1.4349 1.2496 1.6477
## Fracture - Death 1.5717 1.3003 1.8996
## 
## $cancer
##                      HR       L      U
## Well - Fracture  1.0646 0.96748 1.1715
## Well - Death     1.0622 0.96639 1.1675
## Fracture - Death 1.0795 0.94940 1.2273
## 
## $renal
##                      HR       L      U
## Well - Fracture  1.0509 0.73233 1.5080
## Well - Death     1.2291 0.88735 1.7026
## Fracture - Death 1.7631 1.16279 2.6732
## 
## $parkinson
##                      HR       L      U
## Well - Fracture  2.0950 1.36205 3.2224
## Well - Death     1.3059 0.76958 2.2158
## Fracture - Death 1.0982 0.65795 1.8330
## 
## $depression
##                      HR       L      U
## Well - Fracture  1.1124 0.96821 1.2780
## Well - Death     1.0345 0.89889 1.1906
## Fracture - Death 1.0717 0.89250 1.2868
## 
## $`timeperiod[5,Inf)`
##                      HR      L      U
## Well - Fracture  1.3742 1.2683 1.4890
## Well - Death     3.3952 3.0923 3.7278
## Fracture - Death 1.4004 1.1666 1.6810

2.3.4 Transition risks

# Osteopenia:
options(digits = 5)
pmatrix.msm(multi.w2, t = 1, ci = "none", covariates = list(age = 60, Tscore.2 = -1.5), 2)
##             Well Fracture    Death
## Well     0.98867 0.004378 0.006953
## Fracture 0.00000 0.989975 0.010025
## Death    0.00000 0.000000 1.000000
  pmatrix.msm(multi.w2, t = 5, ci = "none", covariates = list(age = 60, Tscore.2 = -1.5), 2)
##             Well Fracture    Death
## Well     0.91057 0.023747 0.065680
## Fracture 0.00000 0.943229 0.056771
## Death    0.00000 0.000000 1.000000
  pmatrix.msm(multi.w2, t = 10, ci = "none", covariates = list(age = 60, Tscore.2 = -1.5), 2)
##             Well Fracture   Death
## Well     0.78473 0.046948 0.16832
## Fracture 0.00000 0.878979 0.12102
## Death    0.00000 0.000000 1.00000
# Osteoporosis:
pmatrix.msm(multi.w2, t = 1, ci = "none", covariates = list(age = 60, Tscore.2 = -2.5), 2)
##             Well  Fracture     Death
## Well     0.99008 0.0028756 0.0070452
## Fracture 0.00000 0.9905687 0.0094313
## Death    0.00000 0.0000000 1.0000000
  pmatrix.msm(multi.w2, t = 5, ci = "none", covariates = list(age = 60, Tscore.2 = -2.5), 2)
##             Well Fracture    Death
## Well     0.91771 0.015672 0.066621
## Fracture 0.00000 0.946515 0.053485
## Death    0.00000 0.000000 1.000000
  pmatrix.msm(multi.w2, t = 10, ci = "none", covariates = list(age = 60, Tscore.2 = -2.5), 2)
##             Well Fracture   Death
## Well     0.79788 0.031186 0.17093
## Fracture 0.00000 0.885750 0.11425
## Death    0.00000 0.000000 1.00000
# Osteoporosis + comorbidities:
pmatrix.msm(multi.w2, t = 1, ci = "none", covariates = list(age = 60, Tscore.2 = -2.5, cvd = 1, diabetes = 1, cancer = 1), 2)
## Warning in factorcov2numeric.msm(covariates, x, mod): Covariates "cvd,
## diabetes" unknown, ignoring
## Warning in msm.fill.pci.covs(x, covariates): Covariate cvd unknown
## Warning in msm.fill.pci.covs(x, covariates): Covariate diabetes unknown
##             Well  Fracture     Death
## Well     0.98946 0.0030593 0.0074823
## Fracture 0.00000 0.9898231 0.0101769
## Death    0.00000 0.0000000 1.0000000
  pmatrix.msm(multi.w2, t = 5, ci = "none", covariates = list(age = 60, Tscore.2 = -2.5, cvd = 1, diabetes = 1, cancer = 1), 2)
## Warning in factorcov2numeric.msm(covariates, x, mod): Covariates "cvd,
## diabetes" unknown, ignoring
## Warning in msm.fill.pci.covs(x, covariates): Covariate cvd unknown
## Warning in msm.fill.pci.covs(x, covariates): Covariate diabetes unknown
##             Well Fracture    Death
## Well     0.91278 0.016611 0.070608
## Fracture 0.00000 0.942390 0.057610
## Death    0.00000 0.000000 1.000000
  pmatrix.msm(multi.w2, t = 10, ci = "none", covariates = list(age = 60, Tscore.2 = -2.5, cvd = 1, diabetes = 1, cancer = 1), 2)
## Warning in factorcov2numeric.msm(covariates, x, mod): Covariates "cvd,
## diabetes" unknown, ignoring
## Warning in msm.fill.pci.covs(x, covariates): Covariate cvd unknown
## Warning in msm.fill.pci.covs(x, covariates): Covariate diabetes unknown
##             Well Fracture   Death
## Well     0.78669 0.032833 0.18048
## Fracture 0.00000 0.877253 0.12275
## Death    0.00000 0.000000 1.00000

2.3.5 Check the model goodness-of-fit

Model validation
prevalence.msm(multi.w2, times = seq(0,20,1), covariates = "mean", ci = "normal")
## $Observed
##    State 1 State 2 State 3 Total
## 0     7269       0       0  7269
## 1     6946     246      74  7266
## 2     6619     428     205  7252
## 3     6345     558     342  7245
## 4     5966     728     529  7223
## 5     5613     866     713  7192
## 6     5242     982     934  7158
## 7     4911    1078    1134  7123
## 8     4573    1154    1359  7086
## 9     4214    1214    1591  7019
## 10    3865    1203    1862  6930
## 11    3499    1226    2108  6833
## 12    3136    1201    2420  6757
## 13    2823    1175    2661  6659
## 14    2468    1122    2953  6543
## 15    2171    1085    3203  6459
## 16    1901    1029    3456  6386
## 17    1603     916    3658  6177
## 18    1163     737    3837  5737
## 19     909     625    4031  5565
## 20     381     299    4153  4833
## 
## $Expected
## $Expected$estimates
##       Well Fracture   Death Total
## 0  7269.00     0.00    0.00  7269
## 1  6942.06   204.70  119.24  7266
## 2  6619.78   392.78  239.44  7252
## 3  6318.54   565.79  360.67  7245
## 4  6018.51   722.95  481.54  7223
## 5  5725.51   864.96  601.54  7192
## 6  5178.83  1035.16  944.01  7158
## 7  4683.58  1176.19 1263.23  7123
## 8  4234.40  1291.21 1560.38  7086
## 9  3811.90  1377.71 1829.39  7019
## 10 3420.39  1439.06 2070.55  6930
## 11 3064.99  1480.38 2287.63  6833
## 12 2754.53  1510.49 2491.98  6757
## 13 2467.05  1522.17 2669.78  6659
## 14 2203.03  1518.02 2821.95  6543
## 15 1976.45  1511.41 2971.14  6459
## 16 1775.92  1499.13 3110.95  6386
## 17 1561.16  1448.03 3167.81  6177
## 18 1317.74  1337.62 3081.64  5737
## 19 1161.68  1285.98 3117.34  5565
## 20  916.88  1103.47 2812.64  4833
## 
## $Expected$ci
## , , 2.5%
## 
##          [,1]    [,2]    [,3] [,4]
##  [1,] 7269.00    0.00    0.00 7269
##  [2,] 6926.03  193.21  110.06 7266
##  [3,] 6589.25  370.66  221.74 7252
##  [4,] 6274.89  534.41  334.19 7245
##  [5,] 5963.13  682.97  446.66 7223
##  [6,] 5659.73  817.27  558.89 7192
##  [7,] 5118.32  991.34  901.15 7158
##  [8,] 4619.61 1129.12 1218.06 7123
##  [9,] 4167.95 1241.66 1511.05 7086
## [10,] 3742.30 1324.17 1775.79 7019
## [11,] 3348.61 1381.66 2011.97 6930
## [12,] 2991.05 1419.01 2224.89 6833
## [13,] 2681.03 1446.22 2423.67 6757
## [14,] 2394.05 1455.94 2597.48 6659
## [15,] 2132.16 1450.19 2746.16 6543
## [16,] 1906.44 1442.54 2891.75 6459
## [17,] 1707.54 1428.61 3028.79 6386
## [18,] 1496.73 1376.51 3087.06 6177
## [19,] 1259.25 1270.36 3005.11 5737
## [20,] 1106.71 1218.98 3040.02 5565
## [21,]  870.76 1043.21 2744.01 4833
## 
## , , 97.5%
## 
##          [,1]    [,2]    [,3] [,4]
##  [1,] 7269.00    0.00    0.00 7269
##  [2,] 6957.27  217.69  128.75 7266
##  [3,] 6648.82  418.45  258.53 7252
##  [4,] 6360.17  602.26  388.80 7245
##  [5,] 6071.44  768.12  518.78 7223
##  [6,] 5788.52  918.39  648.43 7192
##  [7,] 5237.47 1086.11  991.50 7158
##  [8,] 4745.74 1227.08 1313.45 7123
##  [9,] 4299.07 1344.40 1614.12 7086
## [10,] 3876.10 1431.82 1888.27 7019
## [11,] 3487.53 1496.56 2135.14 6930
## [12,] 3132.39 1540.58 2356.79 6833
## [13,] 2825.16 1573.11 2565.09 6757
## [14,] 2537.75 1587.71 2746.06 6659
## [15,] 2272.85 1585.98 2900.24 6543
## [16,] 2045.07 1582.16 3050.91 6459
## [17,] 1843.15 1571.50 3192.64 6386
## [18,] 1624.92 1520.97 3248.50 6177
## [19,] 1375.61 1407.81 3158.33 5737
## [20,] 1216.66 1355.73 3192.70 5565
## [21,]  963.61 1165.26 2878.85 4833
## 
## 
## 
## $`Observed percentages`
##     State 1 State 2 State 3
## 0  100.0000  0.0000  0.0000
## 1   95.5959  3.3856  1.0184
## 2   91.2714  5.9018  2.8268
## 3   87.5776  7.7019  4.7205
## 4   82.5973 10.0789  7.3238
## 5   78.0451 12.0412  9.9138
## 6   73.2327 13.7189 13.0483
## 7   68.9457 15.1341 15.9203
## 8   64.5357 16.2856 19.1787
## 9   60.0370 17.2959 22.6670
## 10  55.7720 17.3593 26.8687
## 11  51.2074 17.9423 30.8503
## 12  46.4111 17.7742 35.8147
## 13  42.3938 17.6453 39.9610
## 14  37.7197 17.1481 45.1322
## 15  33.6120 16.7983 49.5897
## 16  29.7682 16.1134 54.1184
## 17  25.9511 14.8292 59.2197
## 18  20.2719 12.8464 66.8816
## 19  16.3342 11.2309 72.4349
## 20   7.8833  6.1866 85.9301
## 
## $`Expected percentages`
## $`Expected percentages`$estimates
##       Well Fracture   Death
## 0  100.000   0.0000  0.0000
## 1   95.542   2.8173  1.6410
## 2   91.282   5.4161  3.3018
## 3   87.212   7.8093  4.9782
## 4   83.324  10.0090  6.6668
## 5   79.609  12.0266  8.3640
## 6   72.350  14.4616 13.1882
## 7   65.753  16.5126 17.7345
## 8   59.757  18.2220 22.0206
## 9   54.308  19.6283 26.0634
## 10  49.356  20.7656 29.8781
## 11  44.856  21.6651 33.4792
## 12  40.766  22.3545 36.8800
## 13  37.048  22.8589 40.0928
## 14  33.670  23.2007 43.1292
## 15  30.600  23.4001 46.0000
## 16  27.810  23.4752 48.7151
## 17  25.274  23.4422 51.2840
## 18  22.969  23.3156 53.7152
## 19  20.875  23.1084 56.0169
## 20  18.971  22.8320 58.1966
## 
## $`Expected percentages`$ci
## , , 2.5%
## 
##          [,1]    [,2]    [,3]
##  [1,] 100.000  0.0000  0.0000
##  [2,]  95.321  2.6592  1.5147
##  [3,]  90.861  5.1111  3.0576
##  [4,]  86.610  7.3762  4.6126
##  [5,]  82.558  9.4555  6.1838
##  [6,]  78.695 11.3636  7.7710
##  [7,]  71.505 13.8494 12.5894
##  [8,]  64.855 15.8517 17.1004
##  [9,]  58.820 17.5227 21.3244
## [10,]  53.317 18.8655 25.2997
## [11,]  48.320 19.9374 29.0328
## [12,]  43.774 20.7671 32.5610
## [13,]  39.678 21.4033 35.8690
## [14,]  35.952 21.8643 39.0071
## [15,]  32.587 22.1640 41.9709
## [16,]  29.516 22.3338 44.7709
## [17,]  26.739 22.3710 47.4286
## [18,]  24.231 22.2844 49.9767
## [19,]  21.950 22.1433 52.3812
## [20,]  19.887 21.9044 54.6275
## [21,]  18.017 21.5850 56.7764
## 
## , , 97.5%
## 
##          [,1]    [,2]    [,3]
##  [1,] 100.000  0.0000  0.0000
##  [2,]  95.751  2.9960  1.7720
##  [3,]  91.683  5.7702  3.5650
##  [4,]  87.787  8.3128  5.3664
##  [5,]  84.057 10.6343  7.1823
##  [6,]  80.486 12.7696  9.0160
##  [7,]  73.170 15.1734 13.8516
##  [8,]  66.626 17.2271 18.4396
##  [9,]  60.670 18.9726 22.7790
## [10,]  55.223 20.3993 26.9023
## [11,]  50.325 21.5953 30.8100
## [12,]  45.842 22.5461 34.4913
## [13,]  41.811 23.2812 37.9620
## [14,]  38.110 23.8431 41.2383
## [15,]  34.737 24.2394 44.3259
## [16,]  31.662 24.4954 47.2349
## [17,]  28.862 24.6085 49.9944
## [18,]  26.306 24.6231 52.5903
## [19,]  23.978 24.5392 55.0520
## [20,]  21.863 24.3617 57.3711
## [21,]  19.938 24.1105 59.5666
plot.prevalence.msm(multi.w2, mintime = 0, maxtime = 20, legend.pos = c(10, 80), col.obs = "gray", 
                    col.exp = "black")

Diagnostic plots for survival
par(mfrow= c(1,2))
plot.survfit.msm(multi.w2, from = 1, main = "No fracture", range = c(0,20), 
                 ci = "normal", col = "black", col.ci = "black", lty = 2, lwd = 2, 
                 col.surv ="grey", lty.surv = 2, lwd.surv = 1, xlab = "Time (years)")
legend(0.4,0.35, legend = c("Expected", "Expected (95% CI)", "Observed", "Observed (95% CI)"), 
       col = c("black", "black", "grey", "gray"), lty = c(2, 2, 1, 2), lwd = c(2, 1, 2, 1))
plot.survfit.msm(multi.w2, from = 2, main = "Initial fracture", range = c(0,20), 
                 ci = "normal", col = "black", col.ci = "black", lty = 2, lwd = 2, 
                 col.surv ="grey", lty.surv = 2, lwd.surv = 1, xlab = "Time (years)")