Longitudinal Ordinal Analysis for Violet2

Introduction

VIOLET 2 was a randomized clinical trial of seriously ill adults in ICUs to study the effectiveness of vitamin D vs. placebo. It has the advantage of having daily ordinal outcomes assessed for 28 days. The original NEJM paper focused on 1078 patients with confirmed baseline vitamin D deficiency. The analyses presented here use 1352 of the original 1360 randomized patients. Of the 1078 patients, 72 of 539 placebo patients died by day 28, and 92 of 531 vitamin D patients died (proportions of 0.13 and 0.17). Our starting dataset has 1359 patients, with death counts of 89 and 112 (proportions 0.13 and 0.16).

See this and related analyses of the ORCHID study and some additional analysis of VIOLET 2. See this for a study of modeling irregular time spacings with Markov models. Go here for a list of papers about Markov longitudinal modeling of binary and ordinal outcomes.

require(rmsb)
knitrSet(lang='markdown', w=7, h=7, dev='png', fig.path='png/stat-')
alrsimdone <- TRUE
Load(violet2)
d <- violet2
omit <- c('Home', 'Dead', 'In Post Hospital facility')
k <- with(d, sum(day == 1 & bstatus %in% omit))
d <- subset(d, bstatus %nin% omit)
oldstatus <- d$status
levels(d$status) <-
  list(Dead='Dead',
       'Vent/ARDS' = c('On Vent', 'ARDS'),
       'In Hospital/Facility' = c('In Hospital',
                                  'In Post Hospital Facility'),
       Home = 'Home')
table(oldstatus, d$status)
##                            
## oldstatus                    Dead Vent/ARDS In Hospital/Facility  Home
##   Dead                       3749         0                    0     0
##   ARDS                          0       377                    0     0
##   On Vent                       0      2966                    0     0
##   In Hospital                   0         0                 8191     0
##   In Post Hospital Facility     0         0                 2937     0
##   Home                          0         0                    0 19586
d$bstatus        <- droplevels(d$bstatus)
label(d$day)     <- 'Day'
label(d$bstatus) <- 'Baseline Status'
dcf <- d
saveRDS(dcf, 'dcf.rds')
dd <- datadist(d); options(datadist='dd')
options(prType='html')

Descriptive Statistics

7 patients with day 0 status of dead, at home, or in post hospital facility were dropped.

html(describe(d))
d

12 Variables   37856 Observations

id: Reference Id
nmissingdistinct
3785601352
lowest :A01-00208A01-00210A01-00261A01-00337A01-00354
highest:W05-00581W05-00589W05-00590W05-00591W05-00592

bstatus: Baseline Status
image
nmissingdistinct
37828283
 Value             ARDS     On Vent In Hospital
 Frequency         2828       10444       24556
 Proportion       0.075       0.276       0.649
 

day: Day
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
378560280.99914.59.322 2.00 3.00 7.7514.5021.2526.0027.00
lowest : 1 2 3 4 5 , highest: 24 25 26 27 28
status
image
nmissingdistinct
37806504
 Value                      Dead            Vent/ARDS In Hospital/Facility
 Frequency                  3749                 3343                11128
 Proportion                0.099                0.088                0.294
                                
 Value                      Home
 Frequency                 19586
 Proportion                0.518
 

age: SCR age enrolled
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
37856078156.3118.1427.0033.0046.0058.0067.2576.0082.00
lowest : 18 19 20 21 22 , highest: 91 92 93 94 96
treatment
nmissingdistinct
3785602
 Value        Placebo Vitamin D
 Frequency      18676     19180
 Proportion     0.493     0.507
 

lips: Baseline LIPS Score
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
378560280.9965.1753.399 0 0 3 5 7 910
lowest : 0.0 1.0 1.5 2.0 2.5 , highest: 12.0 12.5 13.0 13.5 14.5
charlson: Baseline Charlson Score
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
366801176170.9873.8243.2420023589
lowest : 0 1 2 3 4 , highest: 12 13 14 15 18
 Value          0     1     2     3     4     5     6     7     8     9    10    11
 Frequency   4816  4312  4648  4704  4760  4312  2800  2072  1708   896   644   280
 Proportion 0.131 0.118 0.127 0.128 0.130 0.118 0.076 0.056 0.047 0.024 0.018 0.008
                                         
 Value         12    13    14    15    18
 Frequency    392   168   112    28    28
 Proportion 0.011 0.005 0.003 0.001 0.001
 

sofa: SOFA Total Score
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
378560190.9935.3824.028 0 1 3 5 81012
lowest : 0 1 2 3 4 , highest: 14 15 16 17 18
 Value          0     1     2     3     4     5     6     7     8     9    10    11
 Frequency   2436  3444  3136  3808  4424  3752  3780  2968  2492  2128  1820  1540
 Proportion 0.064 0.091 0.083 0.101 0.117 0.099 0.100 0.078 0.066 0.056 0.048 0.041
                                                     
 Value         12    13    14    15    16    17    18
 Frequency    896   476   308   168   140    84    56
 Proportion 0.024 0.013 0.008 0.004 0.004 0.002 0.001
 

vitD: LCMS Day 0
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
364001456554113.748.279 3.898 5.163 8.04712.60017.80023.31027.300
lowest : 0.0 0.6 0.8 1.3 1.5 , highest: 48.3 56.2 57.0 67.6 68.0
ddeath: Day of Death (29 if alive)
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
378560280.37526.234.918 5132929292929
lowest : 1 2 3 4 5 , highest: 25 26 27 28 29
dhome: Day Discharged to Home (NA if never)
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
2679611060280.9948.1996.851 2 2 4 6111825
lowest : 1 2 3 4 5 , highest: 24 25 26 27 28

Summarize deaths occurring after a patient arrives at home (there were no other events).

u <- subset(d, ! is.na(dhome) & day >= dhome)
np <- length(unique(u$id))
ud <- subset(u, status == 'Dead')
nd <- length(unique(ud$id))
cat(np, 'patients discharged to home\n',
    nd, 'died post discharge\n')
957 patients discharged to home
 29 died post discharge

Event Chart for Random Sample of Patients

A status timeline is shown for each of 65 randomly chosen patients. The patients are sorted so that the worst outcomes at the last follow-up day are at the top. One day is subtract from the time variable to make follow-up start at day 0.

ssamp <- sample(unique(d$id), 65, FALSE)
dr        <- subset(d, id %in% ssamp)
dr        <- subset(dr, day <= ddeath)
dr$id     <- as.integer(as.factor(dr$id))
dr$status <- factor(dr$status, levels=rev(levels(dr$status)))
dr$day    <- dr$day - 1
multEventChart(status ~ day + id, data=dr,
               absorb='Dead', sortbylast = TRUE) +
  theme_classic() +
  theme(legend.position='bottom')

Summary of Successive State Transitions

propsTrans(status ~ day + id, data=d, maxsize=6) +
  theme(legend.position='bottom', axis.text.x=element_text(angle=90, hjust=1))
# For a nice event chart for this type of data see https://livefreeordichotomize.com/2020/05/21/survival-model-detective-1