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)

btxtime

Bayesian Proportional Odds Ordinal Logistic Model

Dirichlet Priors With Concentration Parameter 0.455 for Intercepts

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

Frequencies of Responses

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

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

Compute treatment contrasts (log odds ratios) as a function of time.

k <- contrast(btxtime, list(treatment='Vitamin D', day=1:28),
                       list(treatment='Placebo',   day=1:28))
k

         bstatus day age lips charlson sofa  Contrast      S.E.        Lower
 1   In Hospital   1  58    5        3    5 0.3540591 0.2338975 -0.129725199
 2   In Hospital   2  58    5        3    5 0.3734883 0.2269386 -0.077890340
 3   In Hospital   3  58    5        3    5 0.3928668 0.2221128 -0.066133606
 4*  In Hospital   4  58    5        3    5 0.4119408 0.2195306 -0.030094376
 5*  In Hospital   5  58    5        3    5 0.4304057 0.2190076 -0.006784936
 6*  In Hospital   6  58    5        3    5 0.4479572 0.2200671 -0.002051644
 7*  In Hospital   7  58    5        3    5 0.4642907 0.2220373  0.006342900
 8*  In Hospital   8  58    5        3    5 0.4791018 0.2241803  0.021734510
 9*  In Hospital   9  58    5        3    5 0.4920858 0.2258263  0.022204824
10*  In Hospital  10  58    5        3    5 0.5029385 0.2265063  0.069770881
11   In Hospital  11  58    5        3    5 0.5113711 0.2261697  0.081901485
12*  In Hospital  12  58    5        3    5 0.5171591 0.2254095  0.079327087
13*  In Hospital  13  58    5        3    5 0.5200937 0.2249417  0.092832635
14*  In Hospital  14  58    5        3    5 0.5199662 0.2252829  0.092554319
15*  In Hospital  15  58    5        3    5 0.5165680 0.2265961  0.075674492
16*  In Hospital  16  58    5        3    5 0.5096903 0.2286418  0.079243684
17*  In Hospital  17  58    5        3    5 0.4991245 0.2308418  0.057896195
18*  In Hospital  18  58    5        3    5 0.4846618 0.2324388  0.034690217
19*  In Hospital  19  58    5        3    5 0.4662043 0.2328201 -0.022454378
20*  In Hospital  20  58    5        3    5 0.4440974 0.2320473 -0.048543807
21*  In Hospital  21  58    5        3    5 0.4187972 0.2306030 -0.065062469
22*  In Hospital  22  58    5        3    5 0.3907600 0.2291814 -0.083756952
23*  In Hospital  23  58    5        3    5 0.3604420 0.2285721 -0.093189231
24*  In Hospital  24  58    5        3    5 0.3282992 0.2295390 -0.131587359
25*  In Hospital  25  58    5        3    5 0.2947879 0.2326931 -0.151226537
26*  In Hospital  26  58    5        3    5 0.2603643 0.2383827 -0.212094615
27*  In Hospital  27  58    5        3    5 0.2254845 0.2466329 -0.251975767
28*  In Hospital  28  58    5        3    5 0.1905286 0.2571906 -0.295994535
         Upper Pr(Contrast>0)
 1   0.7925092         0.9318
 2   0.8223911         0.9455
 3   0.8172218         0.9552
 4*  0.8372061         0.9630
 5*  0.8594823         0.9698
 6*  0.8690354         0.9750
 7*  0.8790902         0.9782
 8*  0.9041111         0.9815
 9*  0.9084629         0.9842
10*  0.9594916         0.9860
11   0.9675514         0.9880
12*  0.9614458         0.9898
13*  0.9749873         0.9895
14*  0.9733888         0.9888
15*  0.9639879         0.9870
16*  0.9742668         0.9835
17*  0.9649165         0.9802
18*  0.9499455         0.9778
19*  0.8948266         0.9712
20*  0.8629274         0.9658
21*  0.8381119         0.9600
22*  0.8120379         0.9490
23*  0.7995944         0.9382
24*  0.7616194         0.9175
25*  0.7460600         0.8935
26*  0.7064714         0.8582
27*  0.7005702         0.8185
28*  0.6997217         0.7748

Redundant contrasts are denoted by *

Intervals are 0.95 highest posterior density intervals
Contrast is the posterior mean 

There is a high probability of a treatment effect from day 1 - day 24.

Bayesian Partial Proportional Odds Random Effects Model

Again analyze the “death not carried forward” dataset. Now we relax the proportional hazards assumption with regard to time and treatment. We use the Peterson and Harrell (1990) partial proportional odds model (PPO) which is akin to a Cox model with time-dependent covariates. In general in a longitudinal study one might expect time to not satisfy the PO assumption because the mix of events often changes over time. The departure from PO for time is assumed to be linear whereas the main effect of time uses a 5-knot restricted cubic spline function.

# Check against vgam when clustering is ignored
# Good agreement
if(FALSE) {
  require(VGAM)
  st <- ordered(d$status)
  f <- vgam(st ~ age + treatment, data=d,
            cumulative(reverse=TRUE, parallel=FALSE))
  print(coef(f))
  b <- blrm(status ~ age + treatment, ~ age + treatment, data=d, refresh=100)
  b
}

# Check variance of random effects against mixor
# Don't bother.  Program would not finish even after 2 hours
# require(mixor)
# d$st <- ordered(d$status)
# f <- mixor(st ~ age + day + treatment, id=id, link='logit', data=d)
#f
# sqrt(diag(vcov(f)))
# 

bncppo <- blrm(status ~ bstatus + rcs(age, 4) + rcs(day, 5) +
                        treatment + cluster(id),
                        ~ day + treatment,  # PPO part of model
                                         priorsdppo=0.5, seed=5,
                     data=d, loo=FALSE, iter=3000, refresh=100,
                                         ppairs='bncppo-pairs.png', file='bncppo.rds')

stanDx(bncppo)

Iterations: 3000 on each of 4 chains, with 6000 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                                  480 1.005
y>=In Hospital/Facility                       213 1.021
y>=Home                                       166 1.028
bstatus=On Vent                               226 1.007
bstatus=In Hospital                           256 1.019
age                                           149 1.041
age'                                          225 1.030
age''                                         230 1.017
day                                          3887 1.000
day'                                         6230 1.000
day''                                        7603 1.000
day'''                                       8151 1.000
treatment=Vitamin D                           166 1.024
day:y>=In Hospital/Facility                   367 1.010
treatment=Vitamin D:y>=In Hospital/Facility  8989 1.000
day:y>=Home                                  7463 1.000
treatment=Vitamin D:y>=Home                  9828 1.000
sigmag                                       7349 1.000

stanDxplot(bncppo)

# Tx effects
# Default priors for ppo: -0.1338 0.4077 0.5113 prob 0.35 0.83 0.91
# sd=0.75                 -0.0565 0.3334 0.4256 prob 0.44 0.84 0.91
# sd=0.5                  -0.0109 0.2714 0.3550 prob 0.48 0.85 0.93
# sd=0.5 t prior          -0.0267 0.2783 0.3599 prob 0.46 0.82 0.90 s=4.61
# sd=0.3                   0.0933 0.1746 0.2459 prob 0.63 0.84 0.90
# Above results for t prior no better convergence than with exp

knitr::include_graphics('bncppo-pairs.png')

bncppo

Bayesian Partial Proportional Odds Ordinal Logistic Model

Dirichlet Priors With Concentration Parameter 0.455 for Intercepts

 blrm(formula = status ~ bstatus + rcs(age, 4) + rcs(day, 5) + 
     treatment + cluster(id), ppo = ~day + treatment, data = d, 
     priorsdppo = 0.5, iter = 3000, refresh = 100, loo = FALSE, 
     ppairs = "bncppo-pairs.png", file = "bncppo.rds", seed = 5)
 

Frequencies of Responses

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

      status     bstatus         age         day   treatment cluster(id) 
          50          28           0           0           0           0 
 

Plot posterior density and Compute posterior mean, median, and 0.95 highest posterior density interval for the treatment effect on each category of the response or better.

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

[1] 0.8388333

plot(P(`treatment=Vitamin D`))
P(`treatment=Vitamin D` + `treatment=Vitamin D:y>=In Hospital/Facility` > 0)

[1] 0.8385

plot(P(`treatment=Vitamin D` + `treatment=Vitamin D:y>=In Hospital/Facility`))
P(`treatment=Vitamin D` + `treatment=Vitamin D:y>=Home` > 0)

[1] 0.839

plot(P(`treatment=Vitamin D` + `treatment=Vitamin D:y>=Home`))

Bayes Constrained Partial Proportional Odds Model

The unconstrained partial PO model above allows the treatment effect to be different for every level of the outcome. Here we use the constrained partial PO model to allow treatment to have a special effect on death but to have the same effect on all other endpoints in terms of odds computed on cumulative probabilities. This will allow us to compute the posterior probability that the treatment affects mortality differently than it affects the other outcomes. Unlike above, we do not at the point place skepticism on the differential effect for death. No interaction with time is included in the following model.

The constrained partial PO model’s constraint function cppo never evaluates the departure from PO at the first level of the outcome variable. That is because the probability that Y is the lowest outcome is one minus the probability that it is at or above the second level. So to model a special effect of treatment on death we need to reverse the order of levels in Y.

drev <- transform(d,
                  status = factor(status, rev(levels(status))))
bdcppor <- blrm(status ~ bstatus + rcs(age, 4) + rcs(day, 5) +
                        treatment + cluster(id),
                        ~ day + treatment,  # PPO part of model
                     cppo=function(y) y == 'Dead',
                                         seed=17,
                     data=drev, loo=FALSE, iter=2000, refresh=100,
                                         file='bdcppor.rds')
stanDx(bdcppor)

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>=In Hospital/Facility       97 1.062
y>=Vent/ARDS                 150 1.044
y>=Dead                      282 1.018
bstatus=On Vent               69 1.081
bstatus=In Hospital          103 1.019
age                          127 1.024
age'                         127 1.031
age''                        107 1.038
day                         2140 1.000
day'                        3439 0.999
day''                       4284 1.001
day'''                      4012 1.000
treatment=Vitamin D           90 1.032
day x f(y)                   234 1.009
treatment=Vitamin D x f(y)  2532 1.000
sigmag                      2646 1.000

stanDxplot(bdcppor)

bdcppor

Bayesian Constrained Partial Proportional Odds Ordinal Logistic Model

Dirichlet Priors With Concentration Parameter 0.455 for Intercepts

 blrm(formula = status ~ bstatus + rcs(age, 4) + rcs(day, 5) + 
     treatment + cluster(id), ppo = ~day + treatment, cppo = function(y) y == 
     "Dead", data = drev, iter = 2000, refresh = 100, loo = FALSE, 
     file = "bdcppor.rds", seed = 17)
 

Frequencies of Responses

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

      status     bstatus         age         day   treatment cluster(id) 
          50          28           0           0           0           0 
 

The probability that vitamin D affects death differently than how it affects the non-fatal endpoints is 0.97. Let’s estimate the treatment effect separately for nonfatal and fatal outcomes, in that order.

contrast(bdcppor, list(treatment='Vitamin D'),
                  list(treatment='Placebo'),
         ycut='In Hospital/Facility')

                                  bstatus age  day  Contrast      S.E.
1 status=In Hospital/Facility In Hospital  58 14.5 -0.308613 0.2607144
                                   Lower     Upper Pr(Contrast>0)
1 status=In Hospital/Facility -0.8151987 0.1747911         0.1212

Intervals are 0.95 highest posterior density intervals
Contrast is the posterior mean 

# Note: ycut='Vent/ARDS' would have given the same result as non-PO did
# not kick in until the last outcome level

contrast(bdcppor, list(treatment='Vitamin D'),
                  list(treatment='Placebo'),
         ycut='Dead')

                  bstatus age  day  Contrast      S.E.     Lower     Upper
1 status=Dead In Hospital  58 14.5 0.1194458 0.3425427 -0.558436 0.7623403
              Pr(Contrast>0)
1 status=Dead         0.6368

Intervals are 0.95 highest posterior density intervals
Contrast is the posterior mean 

The probability of benefit of vitamin D on nonfatal outcomes is 1 - 0.12 = 0.88. From the same model the probability of increased mortality is 0.64.

Bayes Binary Longitudinal Model for Death Alone

In what follows we dichotomize the ordinal outcome at death vs. alive and fit a Bayesian longitudinal binary model allowing for a time x treatment interaction, after first fitting a model with no interaction. Because of a lower effective sample size with a binary event, the model is simplified.

d$dead <- 1 * (d$status == 'Dead')
bdeadnoia <- blrm(dead ~ bstatus + rcs(age, 3) + treatment +
                     cluster(id),
                                         seed=13,
                     data=d, loo=FALSE, iter=2000, refresh=100,
                                         file='bdeadnoia.rds')

stanDx(bdeadnoia)

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
Intercept             888 1.005
bstatus=On Vent      2296 1.000
bstatus=In Hospital  2030 1.003
age                  2758 1.001
age'                 3543 1.000
treatment=Vitamin D  3467 1.000
sigmag                651 1.007

stanDxplot(bdeadnoia, 'ALL')

bdeadnoia

Bayesian Logistic Model

Dirichlet Priors With Concentration Parameter 0.541 for Intercepts

 blrm(formula = dead ~ bstatus + rcs(age, 3) + treatment + cluster(id), 
     data = d, iter = 2000, refresh = 100, loo = FALSE, file = "bdeadnoia.rds", 
     seed = 13)
 

        dead     bstatus         age   treatment cluster(id) 
          50          28           0           0           0 
 

The posterior probability that vitamin D increases mortality is 0.836 with this binary logistic model.

Now with the time interaction included.

d$dead <- 1 * (d$status == 'Dead')
bdead <- blrm(dead ~ bstatus + rcs(age, 3) + treatment * rcs(day, 3) +
                     cluster(id),
                                         seed=13,
                     data=d, loo=FALSE, iter=2000, refresh=100,
                                         file='bdead.rds')

stanDx(bdead)

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
Intercept                     21 1.138
bstatus=On Vent              244 1.026
bstatus=In Hospital           71 1.065
age                          204 1.034
age'                        2239 1.000
treatment=Vitamin D          801 1.005
day                           40 1.067
day'                         326 1.021
treatment=Vitamin D * day   1272 1.002
treatment=Vitamin D * day'  1737 1.001
sigmag                        16 1.155

stanDxplot(bdead, 'ALL')

bdead

Bayesian Logistic Model

Dirichlet Priors With Concentration Parameter 0.541 for Intercepts

 blrm(formula = dead ~ bstatus + rcs(age, 3) + treatment * rcs(day, 
     3) + cluster(id), data = d, iter = 2000, refresh = 100, loo = FALSE, 
     file = "bdead.rds", seed = 13)
 

        dead     bstatus         age   treatment         day cluster(id) 
          50          28           0           0           0           0 
 

ggplot(Predict(bdead, day, treatment))

Compute treatment contrasts (log odds ratios) as a function of time.

k <- contrast(bdead, list(treatment='Vitamin D', day=1:28),
                     list(treatment='Placebo',   day=1:28))
k

         bstatus age day     Contrast      S.E.       Lower     Upper
 1   In Hospital  58   1 -0.061995077 0.3118912 -0.63051747 0.5679513
 2   In Hospital  58   2 -0.034768495 0.2781012 -0.55512603 0.5300278
 3*  In Hospital  58   3 -0.007541913 0.2481887 -0.49214380 0.4909035
 4   In Hospital  58   4  0.019673499 0.2237445 -0.39217686 0.4911184
 5*  In Hospital  58   5  0.046821894 0.2066724 -0.33140155 0.4760852
 6*  In Hospital  58   6  0.073836253 0.1982867 -0.29800301 0.4649681
 7*  In Hospital  58   7  0.100649558 0.1985063 -0.28120825 0.4895067
 8*  In Hospital  58   8  0.127194793 0.2056764 -0.27276669 0.5295627
 9*  In Hospital  58   9  0.153404938 0.2172100 -0.27716379 0.5710971
10*  In Hospital  58  10  0.179212976 0.2304379 -0.25741483 0.6293617
11*  In Hospital  58  11  0.204551890 0.2430959 -0.24757546 0.6959855
12*  In Hospital  58  12  0.229354662 0.2534465 -0.24426898 0.7420478
13*  In Hospital  58  13  0.253554273 0.2602500 -0.24087243 0.7671111
14*  In Hospital  58  14  0.277083706 0.2627400 -0.21469462 0.8072200
15*  In Hospital  58  15  0.299897351 0.2608493 -0.18826472 0.8311799
16*  In Hospital  58  16  0.322035234 0.2557943 -0.14843419 0.8565006
17*  In Hospital  58  17  0.343558787 0.2495172 -0.13938872 0.8458842
18*  In Hospital  58  18  0.364529443 0.2444885 -0.08721785 0.8703701
19*  In Hospital  58  19  0.385008635 0.2435442 -0.07064278 0.8740898
20*  In Hospital  58  20  0.405057796 0.2494226 -0.10760083 0.8651527
21*  In Hospital  58  21  0.424738359 0.2640553 -0.08014772 0.9482028
22*  In Hospital  58  22  0.444111757 0.2880380 -0.12805692 0.9852111
23*  In Hospital  58  23  0.463239424 0.3206903 -0.16041240 1.0776938
24*  In Hospital  58  24  0.482182791 0.3605736 -0.16637156 1.2367042
25*  In Hospital  58  25  0.501003292 0.4060202 -0.22888248 1.3578278
26*  In Hospital  58  26  0.519762361 0.4554347 -0.29543608 1.4792000
27*  In Hospital  58  27  0.538511190 0.5074826 -0.38170386 1.5979949
28*  In Hospital  58  28  0.557260020 0.5613802 -0.45143753 1.7464363
     Pr(Contrast>0)
 1           0.4185
 2           0.4440
 3*          0.4815
 4           0.5298
 5*          0.5862
 6*          0.6480
 7*          0.6900
 8*          0.7255
 9*          0.7615
10*          0.7805
11*          0.7972
12*          0.8140
13*          0.8390
14*          0.8578
15*          0.8835
16*          0.9015
17*          0.9200
18*          0.9365
19*          0.9448
20*          0.9500
21*          0.9510
22*          0.9412
23*          0.9275
24*          0.9085
25*          0.8948
26*          0.8718
27*          0.8548
28*          0.8393

Redundant contrasts are denoted by *

Intervals are 0.95 highest posterior density intervals
Contrast is the posterior mean 

GEE Model for Death with Vitamin D Interaction

f <- lrm(dead ~ bstatus + rcs(age, 3) + rcs(day, 3) + rcs(vitD, 3) * treatment,
                     data=d, x=TRUE, y=TRUE)
g <- robcov(f, d$id)
anova(g)

There is mild evidence for an interaction. Look at the contrasts.

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)

Here higher contrasts (higher relative log odds) is bad. Vitamin D treatment seems to move mortality in the wrong direction for patients with low baseline vitamin D, if anything.

Serial Conditional Markov Ordinal Model

Here we use a first order Markov process where the PO model on a given day is conditional on the ordinal outcome on the previous day as well as on the other variables. First we need to create a 1-day lagged outcome status. Unfortunately, since vent and ARDS were pooled for the follow-up status we need to pool them for the baseline status, resulting in some loss of information. Same for in hospital vs. being at another healthcare facility.

d <- dnc
s <- d$bstatus
levels(s) <- list('Vent/ARDS'='ARDS', 'Vent/ARDS'='On Vent', 'In Hospital/Facility'='In Hospital')
table(d$bstatus, s)

             s
              Vent/ARDS In Hospital/Facility
  ARDS             2130                    0
  On Vent          9134                    0
  In Hospital         0                23011

# Now expand baseline levels to include 1 other levels occurring in follow-up
s <- factor(s, setdiff(levels(d$status), 'Dead'))
table(d$bstatus, s)

             s
              Vent/ARDS In Hospital/Facility  Home
  ARDS             2130                    0     0
  On Vent          9134                    0     0
  In Hospital         0                23011     0

d$bstatus <- s
require(data.table)

setDT(d)
d <- d[order(id, day), ]
d[, pstatus := shift(status), by=id]
d[day == 1, pstatus := bstatus]
d$pstatus <- factor(d$pstatus, setdiff(levels(d$pstatus), 'Dead'))
# View(d[, c('id', 'day', 'bstatus', 'status', 'pstatus')])

Fit a fairly saturated model that allows

• transitions to depend on a flexible function of time
• the effect of the previous state to change with time
• the treatment effect to change with time
• a weird effect where the effect of treatment may differ according to the previous outcome state

dd <- datadist(d); options(datadist='dd')
f <- lrm(status ~ rcs(day, 4) * (pstatus + treatment + rcs(age, 3) + lips + charlson + sofa) +
         pstatus * treatment,
         data=d, maxit=30)
f

Logistic Regression Model

 lrm(formula = status ~ rcs(day, 4) * (pstatus + treatment + rcs(age, 
     3) + lips + charlson + sofa) + pstatus * treatment, data = d, 
     maxit = 30)
 

Frequencies of Responses

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

    status       day   pstatus treatment       age      lips  charlson      sofa 
        50         0        49         0         0         0      1045         0 
 

anova(f)

ggplot(Predict(f, day, pstatus))

There was little evidence for a treatment effect at any time, and little evidence that the effect of treatment on state transitions varied with time. There is also little evidence that the treatment effect differed for different previous states.

The graph above shows that being at home is almost an absorbing state, and that the effect of previous state varies with time. Some of that may be due to some odd data: there are 57 patients on ventilator on day 27 but not on day 28.

For now let’s simplify the model by having previous state as the only factor interacting with time, and removing the interactions with treatment. Compare the variances of parameter estimates with cluster sandwich robust estimates.

f <- lrm(status ~ rcs(day, 4) * pstatus + treatment + rcs(age, 3) + lips + charlson + sofa,
         data=d, maxit=30, x=TRUE, y=TRUE)
g <- robcov(f, d$id)
round(diag(vcov(f)) / diag(vcov(g)), 2)

                        y>=Vent/ARDS              y>=In Hospital/Facility 
                                0.86                                 0.92 
                             y>=Home                                  day 
                                0.92                                 1.54 
                                day'                                day'' 
                                1.83                                 1.82 
        pstatus=In Hospital/Facility                         pstatus=Home 
                                0.95                                 0.70 
                 treatment=Vitamin D                                  age 
                                1.07                                 0.95 
                                age'                                 lips 
                                1.08                                 1.11 
                            charlson                                 sofa 
                                1.10                                 1.10 
  day * pstatus=In Hospital/Facility  day' * pstatus=In Hospital/Facility 
                                1.42                                 1.55 
day'' * pstatus=In Hospital/Facility                   day * pstatus=Home 
                                1.58                                 0.70 
                 day' * pstatus=Home                 day'' * pstatus=Home 
                                0.75                                 0.80 

The two variance estimates are similar, giving an indirect indication that the assumption of conditional independence (conditioning on previous state) is reasonable.

Repeat the model truncating the data on day 27 to see how much of the time x previous state interaction resulted from the drop in number of patients on ventilator on day 28.

fm <- lrm(status ~ rcs(day, 4) * pstatus + treatment + rcs(age, 3) + lips + charlson + sofa,
         data=d, maxit=30, subset=day < 28)
fm

Logistic Regression Model

 lrm(formula = status ~ rcs(day, 4) * pstatus + treatment + rcs(age, 
     3) + lips + charlson + sofa, data = d, subset = day < 28, 
     maxit = 30)
 

Frequencies of Responses

                 Dead            Vent/ARDS In Hospital/Facility 
                  188                 3169                10643 
                 Home 
                18088 
 

    status       day   pstatus treatment       age      lips  charlson      sofa 
        48         0        47         0         0         0      1010         0 
 

anova(fm)

ggplot(Predict(fm, day, pstatus))

The interaction effect was cut in half but is still there. However, to put its strength into context, the time by previous status interaction explains only 0.008 of the total amount of outcome variation explained by the status + status x time interaction combination.

We see that being at home is less of an absorbing state for later times, and that the influence of previously being in Vent/ARDS on the current state lessens with time.

# Go back to the dataset that carried death forward, and analyze day 10
f10 <- lrm(status ~ bstatus + treatment, data=dcf, subset=day==10)

priskmarg <- predict(f10, data.frame(treatment='Placebo', bstatus='In Hospital'),
                 type='fitted')['y>=Home']
priskmarg

  y>=Home 
0.6217469 

# Compare with simple estimate not assuming effect of trt & bstatus are additive in log odds
priskmargEmpirical <-
  with(subset(dcf, day == 10 & treatment == 'Placebo' & bstatus=='In Hospital'),
     mean(status == 'Home', na.rm=TRUE))
# Now add covariates
f10c <- lrm(status ~ bstatus + treatment + rcs(age, 3) + lips + charlson + sofa,
            data=dcf, subset=day==10)

w <- subset(dcf, day==10, select=c(age, lips, charlson, sofa))
w <- cbind(w, treatment='Placebo', bstatus='In Hospital')
prisk <- predict(f10c, w, type='fitted')[, 'y>=Home']
mean(prisk, na.rm=TRUE)

[1] 0.5769679

diff <- abs(prisk - priskmargEmpirical)
w$diff <- diff
i <- ! is.na(prisk) & diff < 0.002
sum(i)

[1] 14

setDT(w)
wi <- w[i, ]
wi <- wi[order(diff), ]
wi

    age lips charlson sofa treatment     bstatus         diff
 1:  55  7.0        3    5   Placebo In Hospital 6.198609e-05
 2:  35  8.5        0    8   Placebo In Hospital 9.360643e-05
 3:  31  4.0        6    4   Placebo In Hospital 1.416378e-04
 4:  62  7.0        4    3   Placebo In Hospital 2.850394e-04
 5:  68  5.5        2    5   Placebo In Hospital 3.229592e-04
 6:  26  5.0        0   10   Placebo In Hospital 4.382999e-04
 7:  66  4.0        5    3   Placebo In Hospital 4.643863e-04
 8:  78  3.0        6    0   Placebo In Hospital 8.246585e-04
 9:  29  8.0        2    6   Placebo In Hospital 8.843973e-04
10:  70  1.0        5    4   Placebo In Hospital 1.029117e-03
11:  70  8.0        4    1   Placebo In Hospital 1.511116e-03
12:  56  9.0        1    6   Placebo In Hospital 1.640537e-03
13:  41  6.0        4    5   Placebo In Hospital 1.847777e-03
14:  82  1.5        5    1   Placebo In Hospital 1.907995e-03

# Best match is close to the center
repcov <- wi[1, 1:4]
# For some reason when evaluating the result of Function(fit) with e.g. age=repcov$age
# the arithmetic was wrong.  Need to remove data.frame/data.table and label attributes
repcov <- lapply(repcov, function(x) {x <- unclass(x); attr(x, 'label') <- NULL; x})

# Compute median baseline covariates, using only the first record per patient
co <- all.vars(~ age + lips + charlson + sofa)
setDT(dcf)
m <- dcf[day == 10, sapply(.SD, function(x) list(median = median(x, na.rm=TRUE))),
       .SDcols=co]
m

   age.median lips.median charlson.median sofa.median
1:         58           5               3           5

linpred <- Function(fm)
# Get intercepts
alpha <- coef(fm)[1:3]
# Create a function that computes expit(linpred()) at representative covariate values
# for a given intercept
# This is the conditional probability that Y >= y given previous y = pstatus
pcond <- function(y, pstatus, treatment, day) 
  plogis(alpha[paste0('y>=', y)] + linpred(day=day, pstatus=pstatus, treatment=treatment,
                      age=repcov$age, lips=repcov$lips, charlson=repcov$charlson,
                      sofa=repcov$sofa))
# Test:
pp <- pcond('Home', 'In Hospital/Facility', treatment='Vitamin D', day=1:27)
round(pp, 3)

 [1] 0.080 0.088 0.095 0.103 0.110 0.116 0.121 0.122 0.120 0.114 0.106 0.096
[13] 0.086 0.077 0.068 0.061 0.055 0.051 0.048 0.047 0.046 0.047 0.048 0.050
[25] 0.052 0.054 0.057

# Compare with proportions of same transitions ignoring covariates and treatment
prop <- numeric(27)
for(t in 1:27) {
  i <- with(d, day == t & pstatus == 'In Hospital/Facility')
  cat(sum(is.na(i)), sum(i, na.rm=TRUE), ' ')
  prop[t] <- mean(d$status[i] == 'Home', na.rm=TRUE)
}

1 877  1 877  1 860  1 809  1 727  1 634  1 586  2 521  2 469  2 423  2 405  2 387  2 357  2 334  2 311  2 297  2 280  2 275  2 266  2 258  2 251  2 247  2 241  2 232  2 223  2 217  2 211  

prn(list(prop, pp), fi='/tmp/z')
plot(prop, pp, xlab='Sample Proportion Hosp -> Home',
       ylab='Covariate-Adjusted Model Estimate')
prn(0, fi='/tmp/z')
abline(a=0, b=1, col=gray(.85))

# Sample 50,000 patients for each treatment and create outcomes for each day
N <- 50000
initial.state <- 'In Hospital/Facility'
tlev <- unique(d$treatment)
slev <- levels(d$status)
Y <- array(0L, dim=c(N, 27, 2), dimnames=list(NULL, paste('Day', 1 : 27), tlev))
if(alrsimdone) Y <- readRDS('alrsim.rds') else {
set.seed(1)
Sys.time()
for(i in 1 : N) {   # takes about 28m
  if(i %% 100 == 0) cat(i, file='/tmp/z')
  for(trt in tlev) {
    for(t in 1 : 27) {
      ps <- if(t == 1) initial.state else slev[Y[i, t - 1, trt]]
      if(ps == 'Dead') Y[i, t, trt] <- 1   # absorbing state
      else {
        # Get all 3 cumulative probabilities then 4 cell probabilities
        p1 <- pcond('Vent/ARDS',            ps, treatment=trt, day=t)
        p2 <- pcond('In Hospital/Facility', ps, treatment=trt, day=t)
        p3 <- pcond('Home',                 ps, treatment=trt, day=t)
        cellprobs <- c(1. - p1, p1 - p2, p2 - p3, p3)
        Y[i, t, trt] <- sample(1 : 4, 1, prob=cellprobs)
      }
    }
  }
}
Sys.time()
saveRDS(Y, 'alrsim.rds')
}