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))

 Value             ARDS     On Vent In Hospital
 Frequency         2828       10444       24556
 Proportion       0.075       0.276       0.649
 

 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
 

 Value        Placebo Vitamin D
 Frequency      18676     19180
 Proportion     0.493     0.507
 

 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
 

 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
 

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

Category Proportions Over Time and PO Assumption

The following plot shows how the mix of outcomes changes over time stratified also by treatment.

propsPO(status ~ treatment + day, data=subset(d, day < 15), nrow=1) +
  theme(legend.position='bottom', axis.text.x=element_text(angle=90, hjust=1))

propsPO(status ~ treatment + day, data=subset(d, day > 14), nrow=1) +
  theme(legend.position='bottom', axis.text.x=element_text(angle=90, hjust=1))

Ignoring treatment, compare the observed patterns with what would be expected had the proportional odds assumption been satisfied, fit a PO model containing only day as a predictor, and allow the day effect to be a flexible spline function. Then compute odds ratios against day 1 and apply those odds ratios to the observed proportions from day 1. This is for the original dataset with death carried forward.

f <- lrm(status ~ rcs(day, 5), data=d)
g <- Function(f)   # derive predicted log odds as a function of day
# Make another function that gets the OR vs. day=1
dor <- function(day) exp(g(day) - g(1))
propsPO(status ~ day, odds.ratio=dor, data=d) +
  theme(legend.position='bottom')

Carrying death forward results in a devastating violation of the PO assumption with respect to time. For example, the incidence of death keeps increasing while the incidence of vent/ARDS is decreasing. Now we repeat the process but terminating follow-up at death.

u <- subset(d, day <= ddeath)
f <- lrm(status ~ rcs(day, 5), data=u)
g <- Function(f)   # derive predicted log odds as a function of day
dor <- function(day) exp(g(day) - g(1))
propsPO(status ~ day, odds.ratio=dor, data=u) +
  theme(legend.position='bottom')

The non-PO with respect to time is less problematic when follow-up stops at death but continues for nonfatal events. Note that in the limiting case where death is the only event and the time intervals are small, such truncation results in a binary logistic odds ratio for treatment being almost exactly the hazard ratio in a Cox model analysis when interval risks are small. [However the absolute risk estimates from such a model as a function of time are hard to interpret.].

Note: The lower panel of the above graph is dependent on the choice of reference day. The proportions were estimated from day 1 alone. Instead let’s estimate them from day 4 and repeat.

dor <- function(day) exp(g(day) - g(4))
propsPO(status ~ day, odds.ratio=dor, ref=4, data=u) +
  theme(legend.position='bottom')

Now proportional odds holds with respect to time.

Now do a similar analysis for the second most important predictor: baseline status, still truncating at death. ARDS is taken as the reference group for computing proportions in the follow-up status categories for the purpose of applying the PO assumption.

f <- lrm(status ~ bstatus, data=u)
g <- Function(f)   # derive predicted log odds as a function of bstatus
dor <- function(bstatus) exp(g(bstatus) - g('ARDS'))
propsPO(status ~ bstatus, odds.ratio=dor, ref='ARDS', data=u) 

Proportional odds is reasonably well satisfied with regard to baseline status except that for patients in hospital at baseline, without ARDS or vent, there is an overprediction of the proportion in the vent/ARDS category. Because these summaries are confounded with time, let’s repeat for day 3 only.

u3 <- subset(u, day == 3)
f <- lrm(status ~ bstatus, data=u3)
g <- Function(f)   # derive predicted log odds as a function of bstatus
dor <- function(bstatus) exp(g(bstatus) - g('ARDS'))
propsPO(status ~ bstatus, odds.ratio=dor, ref='ARDS', data=u3) 

See here for assessment of PO for other variables.

GEE-Type Proportional Odds Modeling

This analysis uses a working independence proportional odds model, then the variance-covariance matrix is corrected for within-patient correlation using the robust cluster sandwhich covariance estimator. Covariates are allowed to have nonlinear effects through the use of restricted cubic spline functions with 4 knots. A dot charts shows the relative importance of the predictors. Note that the lrm function in the R rms package frames the model in terms of \(P(Y \geq y | X)\) whereas most software uses \(P(Y \leq y | X)\).

There are several patients excluded in the analysis below due to missing baseline vitamin D levels.

f <- lrm(status ~ bstatus + rcs(day, 5) + rcs(age, 5) + rcs(lips, 5) +
             rcs(charlson, 4) + rcs(sofa, 5) + rcs(vitD, 5) + treatment,
         data=d, tol=1e-13, x=TRUE, y=TRUE)
g <- robcov(f, d$id)
plot(anova(g), trans=sqrt)

We reallocate degrees of freedom for the terms in the model so that the predictors with the most potential from the dot chart are given more parameters. This is because complex relationships that are week will not suffer from underfitting. The likelihood ratio and score \(\chi^2\) statistics in the output below should be ignored because they are not adjusted for intra-cluster correlation.

f <- lrm(status ~ bstatus + rcs(day, 5) + rcs(age, 3) + rcs(lips, 3) +
             rcs(charlson, 4) + rcs(sofa, 4) + rcs(vitD, 3) + treatment,
         data=d, x=TRUE, y=TRUE)
g <- robcov(f, d$id)
g

Logistic Regression Model

 lrm(formula = status ~ bstatus + rcs(day, 5) + rcs(age, 3) + 
     rcs(lips, 3) + rcs(charlson, 4) + rcs(sofa, 4) + rcs(vitD, 
     3) + treatment, data = d, x = TRUE, y = TRUE)
 

Frequencies of Responses

                 Dead            Vent/ARDS In Hospital/Facility 
                 3423                 3043                10467 
                 Home 
                18353 
 

Frequencies of Missing Values Due to Each Variable

anova(g)

Partial effect plots below show the relationship between each predictor and log odds, holding other predictors to median or mode.

ggplot(Predict(g))

A nomogram representation of the model is below. The nomogram can be enhanced to show the estimated probability or outcome better than any given level.

n <- nomogram(g, vnames='names')
plot(n)

d14 <- subset(d, day == 14)
f <- lrm(status ~ bstatus + age + lips +
             charlson + sofa + treatment, data=d14, x=TRUE, y=TRUE)
f

 lrm(formula = status ~ bstatus + age + lips + charlson + sofa + 
     treatment, data = d14, x = TRUE, y = TRUE)
 

                 Dead            Vent/ARDS In Hospital/Facility 
                  142                   77                  303 
                 Home 
                  786 
 

    status   bstatus       age      lips  charlson      sofa treatment 
         2         1         0         0        42         0         0 
 

anova(f)

f <- lrm(status ~ rcs(day, 3) * (bstatus + rcs(age, 3) + rcs(lips, 3) +
             rcs(charlson, 3) + rcs(sofa, 3) + treatment),
                 data=d, tol=1e-13, x=TRUE, y=TRUE)
g <- robcov(f, d$id)
anova(g)

days <- c(1, 8, 15, 22)
ggplot(Predict(g, bstatus, day=days))
ggplot(Predict(g, age, day=days))
ggplot(Predict(g, lips, day=days))
ggplot(Predict(g, charlson, day=days))
ggplot(Predict(g, sofa, day=days))
ggplot(Predict(g, treatment, day=days))
ggplot(Predict(g, day, treatment))

Checking Proportional Odds Assumption

Before getting into a more detailed analysis concerning the proportional odds assumption involving multiple predictors, let’s take a different look at time alone (the approach described earlier is more intuitive). For each day, compute all cumulative proportions, take their logit, and see how the spacing between successive logit cumulative proportions varies across days. Plotting symbols are defined in the code. This may be easier to see in a table of simple cell frequencies that is also shown.

d$yn <- as.integer(d$status)   # numeric Y for easy dichotomization
table(d$status, d$yn)

                      
                           1     2     3     4
  Dead                  3749     0     0     0
  Vent/ARDS                0  3343     0     0
  In Hospital/Facility     0     0 11128     0
  Home                     0     0     0 19586

ql <- function(y) {
  p <- min(0.9999, max(0.0001, mean(y, na.rm=TRUE)))
  qlogis(p)
}
sf <- function(y) c(
  'y>=2'=ql(y >= 2),  # box      (alive)
  'y>=3'=ql(y >= 3),  # circle   (not dead or Vent/ARDS)
  'y>=4'=ql(y == 4))  # triangle (discharged)
s <- summary(yn ~ day, fun=sf, continuous=99, data=d)
plot(s, which=1:3, pch=0:4, xlab='logit Cumulative Prob', main='')
# g <- function(y)
#   apply(y, 2, function(u) qlogis(mean(u)))
# s <- summary(cumcategory(yn) ~ day, continuous=99, data=d)
# plot(s, which=2:3, pch=1:5, xlab='logit Cumulative Prob.', main='')
with(d, print(table(day, status)))

    status
day  Dead Vent/ARDS In Hospital/Facility Home
  1    16       421                  877   37
  2    31       319                  860  141
  3    49       258                  809  235
  4    63       226                  727  335
  5    76       202                  634  439
  6    84       169                  586  512
  7    96       153                  521  580
  8   107       130                  469  644
  9   113       123                  423  691
  10  121       115                  405  709
  11  125       104                  387  734
  12  135        93                  357  765
  13  140        87                  334  789
  14  147        82                  311  810
  15  151        79                  297  823
  16  157        78                  280  835
  17  162        72                  275  841
  18  166        69                  266  849
  19  167        67                  258  858
  20  170        66                  251  863
  21  170        66                  247  867
  22  173        63                  241  873
  23  178        62                  232  878
  24  182        62                  223  883
  25  186        61                  217  886
  26  193        59                  211  887
  27  195        57                  204  894
  28  196         0                  226  928

In the following we ignore intra-patient correlation and use an all-linear (except for time) PO model for purposes of assessing the PO assumption on all covariates including time. The assessment is done by comparing \(\beta\)s with those from the PO model overall fit and those from all dichotomizations of Y. The analysis is done twice. First all patient records are considered. Then days after the day of death are removed. We pay more attention to the more important variables: day, baseline status, SOFA, and Charlson (in descending order of outcome variation explained). The day effect is extremely nonlinear in the PO model, and day is modeled as a quadratic effect here so it has two \(\beta\)s.

for(inclall in c(TRUE, FALSE)) {
  cat(if(inclall) '\nAll records\n\n' else '\nRecords until death\n\n')
  u <- if(inclall) d else subset(d, day <= ddeath)
  f <- lrm(status ~ bstatus + pol(day, 2) + age + lips + charlson + sofa +
             treatment, data=u)
  k <- matrix(NA, nrow=4, ncol=9)
  poc <- coef(f)[-(1:3)]
  colnames(k) <- names(poc)
  rownames(k) <- c('PO', paste0('>=', 2:4))
  k[1, ] <- poc
  for(cut in c(2 : 4)) {
    g <- lrm(yn >= cut ~ bstatus + pol(day, 2) + age + lips + charlson +
               sofa + treatment, data=u)
    k[cut, ] <- coef(g)[-1]
  }
  print(round(k, 3))
}

All records

    bstatus=On Vent bstatus=In Hospital    day  day^2    age   lips charlson
PO            0.464               1.465  0.171 -0.004 -0.009 -0.051   -0.113
>=2           0.582               1.053 -0.183  0.004 -0.004 -0.106   -0.207
>=3           0.410               2.334  0.096 -0.003 -0.001 -0.047   -0.160
>=4           0.364               1.170  0.322 -0.008 -0.014 -0.042   -0.105
      sofa treatment=Vitamin D
PO  -0.089               0.154
>=2 -0.123               0.022
>=3 -0.076               0.046
>=4 -0.098               0.205

Records until death

    bstatus=On Vent bstatus=In Hospital   day  day^2    age   lips charlson
PO            0.326               1.600 0.284 -0.006 -0.015 -0.021   -0.045
>=2           0.490               0.772 0.080 -0.001 -0.006 -0.098   -0.177
>=3           0.271               3.868 0.237 -0.005 -0.003 -0.003   -0.044
>=4           0.193               1.001 0.345 -0.008 -0.018 -0.021   -0.049
      sofa treatment=Vitamin D
PO  -0.075               0.217
>=2 -0.100              -0.035
>=3 -0.034               0.102
>=4 -0.085               0.255

Constancy of the \(\beta\)s vertically in the above tables, and hence equality with the PO \(\beta\)s in the first row, means the PO assumption is satisfied. Not that PO is not strictly satisfied for any dataset we encounter, and even when PO is not satisfied, the PO model can still work better than alternative approaches.

When all records are analyzed (even those with duplicated death status post death), the coefficients are slightly more constant, except for the coefficient of day. When death-carried-forward is not used, the coefficients are slightly less constant except for the effect of day. To help explore that, here is an analysis of when deaths are occuring (patients not dying are removed):

w <- subset(d, day == 1 & ddeath < 29)
hist(w$ddeath, breaks=0:28, main='', xlab='Day of Death')

GEE-Like Analysis Repeated Without Carrying Death Forward

The last analysis with cluster correlation correction is repeated after deletion of patient records beyond the day of death. In other words, the death status is not repeated up to day 28 but the patient data are truncated at death.

d <- dnc <- subset(d, day <= ddeath)
f <- lrm(status ~ bstatus + rcs(day, 5) + rcs(age, 3) + rcs(lips, 3) +
             rcs(charlson, 4) + rcs(sofa, 4) + rcs(vitD, 3) + treatment,
         data=d, x=TRUE, y=TRUE)
g <- robcov(f, d$id)
g

Logistic Regression Model

 lrm(formula = status ~ bstatus + rcs(day, 5) + rcs(age, 3) + 
     rcs(lips, 3) + rcs(charlson, 4) + rcs(sofa, 4) + rcs(vitD, 
     3) + treatment, data = d, x = TRUE, y = TRUE)
 

Frequencies of Responses

                 Dead            Vent/ARDS In Hospital/Facility 
                  180                 3043                10467 
                 Home 
                18353 
 

Frequencies of Missing Values Due to Each Variable

anova(g)

ggplot(Predict(g))

Let’s also look at the proportional odds assumption using these truncated records (death not carried forward).

s <- summary(yn ~ day, fun=sf, continuous=99, data=d)
plot(s, which=1:3, pch=0:4, xlab='logit Cumulative Prob', main='')

When the absorbing state is not carried forward, the problem with non-proportional odds almost vanishes.

Finally, look at interactions with time in this model.

f <- lrm(status ~ rcs(day, 3) * (bstatus + rcs(age, 3) + rcs(lips, 3) +
             rcs(charlson, 3) + rcs(sofa, 3) + treatment),
                 data=d, tol=1e-13, x=TRUE, y=TRUE)
g <- robcov(f, d$id)
anova(g)

days <- c(1, 8, 15, 22)
ggplot(Predict(g, bstatus, day=days))
ggplot(Predict(g, age, day=days))
ggplot(Predict(g, lips, day=days))
ggplot(Predict(g, charlson, day=days))
ggplot(Predict(g, sofa, day=days))
ggplot(Predict(g, treatment, day=days))
ggplot(Predict(g, day, treatment))

f <- lrm(status ~ bstatus + rcs(day, 5) + rcs(age, 3) + rcs(lips, 3) +
             rcs(charlson, 4) + rcs(sofa, 4) + rcs(vitD, 4) * treatment,
         data=d, x=TRUE, y=TRUE)
g <- robcov(f, d$id)
g

 lrm(formula = status ~ bstatus + rcs(day, 5) + rcs(age, 3) + 
     rcs(lips, 3) + rcs(charlson, 4) + rcs(sofa, 4) + rcs(vitD, 
     4) * treatment, data = d, x = TRUE, y = TRUE)
 

                 Dead            Vent/ARDS In Hospital/Facility 
                  180                 3043                10467 
                 Home 
                18353 
 

anova(g)

ggplot(Predict(g, vitD, treatment))

k <- contrast(g, list(treatment='Vitamin D', vitD=0:65),
                 list(treatment='Placebo', vitD=0:65))
k <- as.data.frame(k[c('vitD', 'Contrast', 'Lower', 'Upper')])
ggplot(k, aes(x=vitD, y=Contrast)) + geom_point() + 
  geom_errorbar(aes(ymin=Lower, ymax=Upper), width=0)

Bayesian Random Effects Model Without Carrying Death Forward

Repeat the last model but using a Bayesian approach with normally distributed random effects for patient level intercepts. This is a compound symmetry type of correlation structure within patient. Before proceeding witih the full model try to find the smallest model that exhibits convergence problems that will be seen later with the full model.

stanSet()

# This is fine:
bnocsnoc <- blrm(status ~ age, data=d, loo=FALSE, refresh=100,
                 file='bnocsnoc.rds')
stanDx(bnocsnoc)

Iterations: 2000 on each of 4 chains, with 4000 posterior distribution samples saved

For each parameter, n_eff is a crude measure of effective sample size
and Rhat is the potential scale reduction factor on split chains
(at convergence, Rhat=1)

                        n_eff  Rhat
y>=Vent/ARDS             4146 1.000
y>=In Hospital/Facility  4713 1.000
y>=Home                  4700 0.999
age                      4419 0.999

stanDxplot(bnocsnoc, which='ALL')

# Problem with this one if rsdmean=0.3 or 0.08, minor prob if 0.005
bnocs <- blrm(status ~ age + cluster(id), psigma=1, rsdsd=0.075,
              data=d, loo=FALSE, refresh=100, file='bnocs.rds')
# Previous: psigma=2 rsdmean=0.005; psigma=1 rsdmean=0.08
stanDx(bnocs)

Iterations: 2000 on each of 4 chains, with 4000 posterior distribution samples saved

For each parameter, n_eff is a crude measure of effective sample size
and Rhat is the potential scale reduction factor on split chains
(at convergence, Rhat=1)

                        n_eff  Rhat
y>=Vent/ARDS              421 1.003
y>=In Hospital/Facility   194 1.011
y>=Home                   173 1.011
age                       199 1.009
sigmag                    385 1.008

stanDxplot(bnocs, which='ALL')
# Problematic even if binary for Home, but got slightly better with
# t(4, 0, 1) prior for sigma SD
# adding conc=0.5 rsdmean made things significanty worse
# b <- blrm(status ~ age + cluster(id), data=d, iter=1000, refresh=100)
# stanDx(b); stanDxplot(b, which='ALL')

bnocs

Bayesian Proportional Odds Ordinal Logistic Model

Dirichlet Priors With Concentration Parameter 0.455 for Intercepts

 blrm(formula = status ~ age + cluster(id), data = d, psigma = 1, 
     rsdsd = 0.075, refresh = 100, loo = FALSE, file = "bnocs.rds")
 

Frequencies of Responses

                 Dead            Vent/ARDS In Hospital/Facility 
                  196                 3343                11128 
                 Home 
                19586 
 

      status         age cluster(id) 
          50           0           0 
 

bnocarry <- blrm(status ~ bstatus + rcs(day, 5) + rcs(age, 3) +
                 rcs(lips, 3) + rcs(charlson, 4) + rcs(sofa, 4) + treatment +
                         cluster(id), data=d, loo=FALSE, refresh=50, file='bnocarry.rds')
bnocarry

Bayesian Proportional Odds Ordinal Logistic Model

Dirichlet Priors With Concentration Parameter 0.455 for Intercepts

 blrm(formula = status ~ bstatus + rcs(day, 5) + rcs(age, 3) + 
     rcs(lips, 3) + rcs(charlson, 4) + rcs(sofa, 4) + treatment + 
     cluster(id), data = d, refresh = 50, loo = FALSE, file = "bnocarry.rds")
 

Frequencies of Responses

                 Dead            Vent/ARDS In Hospital/Facility 
                  189                 3169                10862 
                 Home 
                18988 
 

      status     bstatus         day         age        lips    charlson 
          50          28           0           0           0        1045 
        sofa   treatment cluster(id) 
           0           0           0 
 

a <- anova(bnocarry)
a

stanDx(bnocarry)

Iterations: 2000 on each of 4 chains, with 4000 posterior distribution samples saved

For each parameter, n_eff is a crude measure of effective sample size
and Rhat is the potential scale reduction factor on split chains
(at convergence, Rhat=1)

                        n_eff  Rhat
y>=Vent/ARDS              336 1.008
y>=In Hospital/Facility   138 1.026
y>=Home                   101 1.038
bstatus=On Vent           110 1.042
bstatus=In Hospital       159 1.010
day                      1357 1.000
day'                     2714 1.001
day''                    3691 1.001
day'''                   4264 1.000
age                       141 1.032
age'                      163 1.022
lips                      121 1.038
lips'                      52 1.076
charlson                  126 1.027
charlson'                 141 1.015
charlson''                122 1.028
sofa                      143 1.024
sofa'                     140 1.024
sofa''                    159 1.017
treatment=Vitamin D       130 1.039
sigmag                    264 1.005

stanDxplot(bnocarry)
ggplot(Predict(bnocarry))

plot(a)

Compute the posterior probability that treatment raises the chance of a patient having a given outcome or better (this is already in the coefficient table above).

P <- PostF(bnocarry, name='orig')
P(`treatment=Vitamin D` > 0)

[1] 0.9575

For comparison, fit a similar Bayesian model adjusting only for time, and one not adjusting for any covariates.

bnocarry2 <- blrm(status ~ rcs(day, 5) + treatment +
                  cluster(id), data=d, loo=FALSE, refresh=50, file='bnocarry2.rds')
bnocarry2

Bayesian Proportional Odds Ordinal Logistic Model

Dirichlet Priors With Concentration Parameter 0.455 for Intercepts

 blrm(formula = status ~ rcs(day, 5) + treatment + cluster(id), 
     data = d, refresh = 50, loo = FALSE, file = "bnocarry2.rds")
 

Frequencies of Responses

                 Dead            Vent/ARDS In Hospital/Facility 
                  196                 3343                11128 
                 Home 
                19586 
 

      status         day   treatment cluster(id) 
          50           0           0           0 
 

bnocarry3 <- blrm(status ~ treatment + cluster(id),
                  data=d, loo=FALSE, refresh=50, file='bnocarry3.rds')
bnocarry3

Bayesian Proportional Odds Ordinal Logistic Model

Dirichlet Priors With Concentration Parameter 0.455 for Intercepts

 blrm(formula = status ~ treatment + cluster(id), data = d, refresh = 50, 
     loo = FALSE, file = "bnocarry3.rds")
 

Frequencies of Responses

                 Dead            Vent/ARDS In Hospital/Facility 
                  196                 3343                11128 
                 Home 
                19586 
 

      status   treatment cluster(id) 
          50           0           0 
 

Adjusting only for time and not other baseline covariates, the between-patient outcome heterogeneity is greater as makes sense, and the treatment effect remains diminished. When not even adjusting for time, the value of \(\sigma_\gamma\) ironically becomes lower than from the model with full covariate adjustment, and the treatment effect becomes very small.

Bayesian PO Random Effects Model With Time by Treatment Interaction

We go back to adjusting for the full set of covariates, and allow time to nonlinearly interact with treatment.

btxtime <- blrm(status ~ bstatus + treatment * rcs(day, 4) + rcs(age, 3) +
                 rcs(lips, 3) + rcs(charlson, 4) + rcs(sofa, 4) +
                         cluster(id), data=d, loo=FALSE, refresh=50, file='btxtime.rds')
stanDx(btxtime)

Iterations: 2000 on each of 4 chains, with 4000 posterior distribution samples saved

For each parameter, n_eff is a crude measure of effective sample size
and Rhat is the potential scale reduction factor on split chains
(at convergence, Rhat=1)

                            n_eff  Rhat
y>=Vent/ARDS                  309 1.013
y>=In Hospital/Facility       149 1.022
y>=Home                       120 1.023
bstatus=On Vent               142 1.027
bstatus=In Hospital           143 1.020
treatment=Vitamin D           154 1.011
day                          1979 1.005
day'                         3281 1.000
day''                        3943 1.001
age                           131 1.026
age'                          119 1.031
lips                           67 1.064
lips'                         149 1.012
charlson                      144 1.038
charlson'                     155 1.033
charlson''                     79 1.066
sofa                          168 1.017
sofa'                         164 1.011
sofa''                        104 1.053
treatment=Vitamin D * day    3919 1.000
treatment=Vitamin D * day'   3880 1.000
treatment=Vitamin D * day''  4018 1.000
sigmag                        298 1.010

stanDxplot(btxtime)