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$id     <- as.integer(as.factor(dr$id))
dr$status <- factor(dr$y, levels=7:1)
dr$day    <- dr$day - 1
multEventChart(status ~ day + id, data=dr,
               absorb='1', sortbylast = TRUE) +
  theme_classic() +
  theme(legend.position='bottom')

Summary of Successive State Transitions

gg <- propsTrans(y ~ day + id, data=dall, maxsize=6, labels=ordLevels) +
        theme(legend.position='bottom',
        axis.text.x=element_text(angle=90, hjust=1))
ggp(gg)

# 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

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

gg <- propsPO(y ~ trt + day, nrow=1, data=dall) +
        theme(legend.position='bottom',
                      axis.text.x=element_text(angle=90, hjust=1))
ggp(gg)

Proportional Odds Assumption

The treatment is too weak to detect any non-proportional odds. Let’s check the PO assumption on the two strongest variables: gap and yprev. We omit important interactions between them and make day linear to be able to limit the number of regression coefficient estimates to examine as a function of y-cutoff. Expect results for the first cutoff, which corresponds to predicting deaths, to be unstable due to there being only 50 deaths.

n <- c('yprev=3', 'yprev=4', 'yprev=5', 'yprev=6', 'yprev=7', 'day', 'gap')
z <- se <- matrix(NA, nrow=6, ncol=length(n))
dimnames(z) <- dimnames(se) <- list(paste0('y>=', 2:7), n)
for(ycut in 2 : 7) {
    g <- lrm(y >= ycut ~ gap + yprev + day + trt + rcs(age, 3) + race +
                     sofa_nogcs + sex,
                     data=d, tol=1e-14, maxit=30)
    z [ycut - 1, n] <- coef(g)[n]
    se[ycut - 1, n] <- sqrt(diag(vcov(g)[n, n])) 
}
round(z, 2)

     yprev=3 yprev=4 yprev=5 yprev=6 yprev=7   day   gap
y>=2    0.41    1.70    2.76    2.71   10.61 -0.30  0.35
y>=3    3.29    5.09    6.43    6.20   13.45 -0.05  0.11
y>=4    1.20    4.82    6.59    6.72   12.62  0.24 -0.32
y>=5    1.94    3.35    6.70    8.10    8.02  0.41 -0.53
y>=6    2.58    3.82    5.01    8.42    7.42  0.64 -0.93
y>=7    2.32    3.06    3.75    4.50    6.08  0.28 -0.35

round(se, 2)

     yprev=3 yprev=4 yprev=5 yprev=6 yprev=7  day  gap
y>=2    0.47    0.45    0.69    0.64   27.72 0.24 0.47
y>=3    0.24    0.28    0.54    0.54   30.24 0.10 0.20
y>=4    0.26    0.27    0.41    0.55   20.70 0.08 0.16
y>=5    0.42    0.39    0.42    0.54    1.07 0.07 0.13
y>=6    0.55    0.51    0.52    0.56    0.64 0.07 0.14
y>=7    1.08    1.02    1.02    1.02    1.03 0.11 0.20

The effect of day and gap are roughly linear in the y-cutoff.
yprev may be predicting death differently from how it predicts the nonfatal outcomes.

The VGAM package vglm function would not converge for a partial PO model that allowed for non-PO in yprev.

Computing Unconditional Estimates

One can uncondition the previous model to obtain unconditional treatment effects, although these will be covariate- and time-dependent. We specify covariate settings (here we use representative values for each covariate, just for convenience, using a method described below) and vary the treatment assignment. Note the following.

1. The initial state (baseline ordinal status) is always conditioned on, so is not marginalized out. Calculations can be done separately for each initial state.
2. Given this initial state, use the fitted model to get the conditional distribution of Y at (post-randomization) day 1 given the initial state, covariate settings, and treatment.
3. Likewise for later days.
4. Death is an absorbing state so once death occurs all future status values are also set to death.
5. Were time to be continuous one can find the first time for which the probability of being in state y or worse exceeds 0.5, which would define the median time until reaching that state. For the ORCHID study, gaps between measurements are too long for this unless time-interpolation is used. It may be better to focus on estimating the probability of being in state y or worse as a function of days since randomization, treatment, and covariates (covariates include the initial state).
6. This approach does not factor in uncertainties in parameter estimates, so provides only a point estimate. The simple Bayesian solution to this is provided later.

Start with computation of a set of representative covariate values that will be used when getting estimates related to transformed parameters. This is done by using the global transition model fitted above, evaluating it at day 14 for a previous state of y=4, and finding covariate settings that give predicted values that are close to the unadjusted estimates for placebo. Of the qualifying settings, one that is subjectively closer to the center of the data distribution is chosen. For simplicity we drop weak categorical predictors sex and race. To better fit we allow time to have a nonlinear effect. First let’s see the distribution of day by yprev.

ggplot(d, aes(x=day)) + geom_histogram() + facet_wrap(~ yprev)

# For yprev=7 there are only 2 distinct assessment days (14, 28)
# So need interaction to be restricted to be linear
f <- lrm(y ~ gap * yprev + rcs(day, 4) + trt + rcs(age, 3) + sofa_nogcs,
         data=d)
newdat <- subset(d, day == 1)
newdat <- transform(newdat, day=14, yprev='4', gap=1)
xb1 <- predict(f, newdat)  # uses first intercept by default

# Do the same when covariates are omitted
fu <- lrm(y ~ gap * yprev + rcs(day, 4), data=d)
xb2 <- predict(fu, newdat[1, ])
diff <- abs(xb1 - xb2)
w <- newdat[, c('age', 'sofa_nogcs', 'trt')]
w$diff <- abs(xb1 - xb2)
setDT(w)
wi <- w[order(diff), ]
# Compute mean and median baseline covariates
co <- all.vars(~ age + sofa_nogcs)
setDT(newdat)
h <- function(x) c(mean=mean(x), median=median(x))
u <- function(x) { s <- sapply(x, h); w <- as.list(s)
                   names(w) <- outer(rownames(s), colnames(s), paste, sep='.'); w}
newdat[, u(.SD), .SDcols=co]

   mean.age median.age mean.sofa_nogcs median.sofa_nogcs
1: 56.85386         57        2.814196                 2

# Seek age=57, sofa 2 or 3
wi[1:30, ]

    age sofa_nogcs                trt        diff
 1:  47          7 hydroxychloroquine 0.001365465
 2:  63          5 hydroxychloroquine 0.001857581
 3:  76          3 hydroxychloroquine 0.002198086
 4:  75          3            placebo 0.002726674
 5:  88          1 hydroxychloroquine 0.005283669
 6:  70          4 hydroxychloroquine 0.006873129
 7:  82          2            placebo 0.007587458
 8:  48          7 hydroxychloroquine 0.008507628
 9:  63          5            placebo 0.010995304
10:  45          7            placebo 0.011599411
11:  75          3 hydroxychloroquine 0.011864397
12:  75          3 hydroxychloroquine 0.011864397
13:  83          2 hydroxychloroquine 0.012547716
14:  56          6            placebo 0.014215440
15:  64          5 hydroxychloroquine 0.014776646
16:  64          5 hydroxychloroquine 0.014776646
17:  77          3 hydroxychloroquine 0.016282989
18:  77          3 hydroxychloroquine 0.016282989
19:  48          7            placebo 0.017645351
20:  89          1            placebo 0.017952035
21:  89          1            placebo 0.017952035
22:  89          1            placebo 0.017952035
23:  87          1 hydroxychloroquine 0.019381651
24:  68          4 hydroxychloroquine 0.020411089
25:  71          4 hydroxychloroquine 0.020655362
26:  61          5 hydroxychloroquine 0.023464890
27:  61          5 hydroxychloroquine 0.023464890
28:  80          2 hydroxychloroquine 0.029746227
29:  67          4 hydroxychloroquine 0.033890653
30:  72          4 hydroxychloroquine 0.034516065
    age sofa_nogcs                trt        diff

# There are no good compromises so use marginal representative values
repcov <- list(age=57, sofa_nogcs=3)

There were no representative covariate combinations that were sufficiently close to unadjusted estimates and close to the center of covariate distributions, so we chose marginally central values.

Compute \(\Pr(Y_s \geq y | X, Y_0)\) analytically for a single value of \((y, Y_0, X, s)\), using this recursive formula where \(k\) is the number of outcome states and \(t\) is time:

\[\Pr(Y_{t} = y | X) = \sum_{j=1}^{k} \Pr(Y_{t} = y | X, Y_{t-1} = j) \Pr(Y_{t-1} = j | X)\] Note that recent updates to the R Hmisc package make the uprob function below obsolete. The soprobMarkovOrd function will compute state occupancy probabilities for given covariate settings. For Bayesian models Hmisc has soprobMarkovOrdm.

# Compute unconditional probabilities
# Output is a matrix with s rows corresponding to t=times(1), times(2), ...
# Each row has columns corresponding to the probabilities of being in each of the states
# times(1) must be 1
# 
# fit is the model fit which must contain the variable named day and the
# previous state, which must be named by the value of pvarname
# data is a data frame or data table that contains covariate settings
# other than day and what's in pvarname.  data has one row.
# times is the vector of days for which to evaluate state occupancy probabilities
# ylevels are the ordered levels of the outcome variable
# Set gap=TRUE if predictor gap is in the model and gap needs to be computed
# At t=times[1] gap is taken as times[1]

uprob <- function(fit, data, times, ylevels, pvarname='yprev', gap=FALSE) {
  k <- length(ylevels)
  s <- length(times)
  P <- matrix(NA, nrow=s, ncol=k)
  colnames(P) <- as.character(ylevels)
  rownames(P) <- as.character(times)
  # Never uncondition on initial state
  data$day <- times[1]
  if(gap) data$gap <- times[1]
  p <- predict(fit, data, type='fitted.ind')
  P[1, ] <- p
  # Matrix of conditional probabilities of Y conditioning on previous time Y
  # Columns = k conditional probabilities conditional on a single previous state
  # Rows    = all possible previous states
  cp <- matrix(NA, nrow=k, ncol=k, 
               dimnames=list(paste('t-1', ylevels), paste('t', ylevels)))
  data <- as.list(data)
  data[[pvarname]] <- as.character(ylevels[-1])
  # Above: don't request estimates after absorbing state
  edata <- expand.grid(data)
  for(it in 2 : s) {
    edata$day <- times[it]
    if(gap) edata$gap <- times[it] - times[it - 1]
    # edata$gap <- c(1,1,1,1,3,3,3,7,7,7,7,7,7,7,rep(14, 14))[times[it]]
    # y=1 = death = absorbing state
    cp <- rbind(c(1., rep(0., k - 1)), predict(fit, edata, type='fitted.ind'))
    # Compute unconditional probabilities of being in all possible states
    # at current time t
    P[it, ] <- t(cp) %*% P[it - 1, ]
    }
  P
}

dat <- data.frame(yprev='4', age=repcov$age, sofa_nogcs=repcov$sofa_nogcs)
  
times <- c(1, 2, 3, 4, 7, 14, 28)
tlev  <- levels(d$trt)

U <- NULL
for(trt in tlev) {   # separately for treatments
    dat$trt <- trt
    u <- uprob(f, dat, times, 1:7, gap=TRUE)
  cat('\n', trt, '\n\n')
  print(round(u, 3))
  U <- rbind(U, data.frame(trt, day=times, p=u, phome=u[, 6] + u[, 7]))
}

 placebo 

       1     2     3     4     5     6     7
1  0.001 0.037 0.172 0.646 0.118 0.024 0.001
2  0.007 0.085 0.157 0.491 0.173 0.077 0.008
3  0.016 0.100 0.130 0.395 0.194 0.141 0.024
4  0.028 0.107 0.117 0.339 0.189 0.184 0.036
7  0.033 0.076 0.077 0.223 0.209 0.282 0.100
14 0.038 0.051 0.038 0.104 0.150 0.394 0.225
28 0.042 0.035 0.014 0.036 0.056 0.374 0.444

 hydroxychloroquine 

       1     2     3     4     5     6     7
1  0.001 0.036 0.171 0.646 0.119 0.025 0.001
2  0.007 0.084 0.156 0.491 0.174 0.078 0.009
3  0.016 0.098 0.129 0.394 0.195 0.142 0.024
4  0.027 0.106 0.116 0.338 0.190 0.186 0.036
7  0.033 0.075 0.076 0.222 0.209 0.284 0.101
14 0.037 0.050 0.038 0.103 0.149 0.395 0.229
28 0.041 0.034 0.013 0.035 0.055 0.373 0.448

Rorchid <- U

Plot the unconditional probabilities of being at home at each time, for representative age and sofa score.

xl <- xlab('Day'); yl <- ylab('Probability of Being at Home')
U$txt <- with(U, paste0(trt, '<br>Day:', day, '<br>', round(phome, 3)))
gg <- ggplot(U, aes(x=day, y=phome, color=trt, label=txt)) +
      geom_line() +
      xl + yl + labs(color='')
ggp(gg)

Repeat but estimate probabilities at days not observed in the data.

U$type <- 'Actual Times'
for(trt in tlev) {   # separately for treatments
    dat$trt <- trt
    u <- uprob(f, dat, 1:28, 1:7, gap=TRUE)
  cat('\n', trt, '\n\n')
  print(round(u, 3))
  lab <- paste(trt, 'day 1-28')
  U <- rbind(U, data.frame(trt, day=1:28, p=u, phome=u[, 6] + u[, 7],
                                                 type='Days 1-28', txt=''))
}

 placebo 

       1     2     3     4     5     6     7
1  0.001 0.037 0.172 0.646 0.118 0.024 0.001
2  0.007 0.085 0.157 0.491 0.173 0.077 0.008
3  0.016 0.100 0.130 0.395 0.194 0.141 0.024
4  0.028 0.107 0.117 0.339 0.189 0.184 0.036
5  0.046 0.126 0.119 0.317 0.179 0.184 0.030
6  0.084 0.162 0.138 0.309 0.162 0.129 0.014
7  0.170 0.216 0.168 0.275 0.111 0.056 0.004
8  0.342 0.264 0.164 0.173 0.042 0.013 0.001
9  0.609 0.231 0.091 0.058 0.008 0.002 0.000
10 0.867 0.102 0.022 0.009 0.001 0.000 0.000
11 0.980 0.017 0.002 0.001 0.000 0.000 0.000
12 0.999 0.001 0.000 0.000 0.000 0.000 0.000
13 1.000 0.000 0.000 0.000 0.000 0.000 0.000
14 1.000 0.000 0.000 0.000 0.000 0.000 0.000
15 1.000 0.000 0.000 0.000 0.000 0.000 0.000
16 1.000 0.000 0.000 0.000 0.000 0.000 0.000
17 1.000 0.000 0.000 0.000 0.000 0.000 0.000
18 1.000 0.000 0.000 0.000 0.000 0.000 0.000
19 1.000 0.000 0.000 0.000 0.000 0.000 0.000
20 1.000 0.000 0.000 0.000 0.000 0.000 0.000
21 1.000 0.000 0.000 0.000 0.000 0.000 0.000
22 1.000 0.000 0.000 0.000 0.000 0.000 0.000
23 1.000 0.000 0.000 0.000 0.000 0.000 0.000
24 1.000 0.000 0.000 0.000 0.000 0.000 0.000
25 1.000 0.000 0.000 0.000 0.000 0.000 0.000
26 1.000 0.000 0.000 0.000 0.000 0.000 0.000
27 1.000 0.000 0.000 0.000 0.000 0.000 0.000
28 1.000 0.000 0.000 0.000 0.000 0.000 0.000

 hydroxychloroquine 

       1     2     3     4     5     6     7
1  0.001 0.036 0.171 0.646 0.119 0.025 0.001
2  0.007 0.084 0.156 0.491 0.174 0.078 0.009
3  0.016 0.098 0.129 0.394 0.195 0.142 0.024
4  0.027 0.106 0.116 0.338 0.190 0.186 0.036
5  0.045 0.124 0.118 0.316 0.180 0.186 0.030
6  0.083 0.161 0.137 0.309 0.164 0.132 0.015
7  0.167 0.214 0.168 0.277 0.113 0.057 0.005
8  0.337 0.264 0.165 0.176 0.043 0.014 0.001
9  0.604 0.233 0.093 0.060 0.008 0.002 0.000
10 0.863 0.104 0.023 0.009 0.001 0.000 0.000
11 0.979 0.018 0.002 0.001 0.000 0.000 0.000
12 0.999 0.001 0.000 0.000 0.000 0.000 0.000
13 1.000 0.000 0.000 0.000 0.000 0.000 0.000
14 1.000 0.000 0.000 0.000 0.000 0.000 0.000
15 1.000 0.000 0.000 0.000 0.000 0.000 0.000
16 1.000 0.000 0.000 0.000 0.000 0.000 0.000
17 1.000 0.000 0.000 0.000 0.000 0.000 0.000
18 1.000 0.000 0.000 0.000 0.000 0.000 0.000
19 1.000 0.000 0.000 0.000 0.000 0.000 0.000
20 1.000 0.000 0.000 0.000 0.000 0.000 0.000
21 1.000 0.000 0.000 0.000 0.000 0.000 0.000
22 1.000 0.000 0.000 0.000 0.000 0.000 0.000
23 1.000 0.000 0.000 0.000 0.000 0.000 0.000
24 1.000 0.000 0.000 0.000 0.000 0.000 0.000
25 1.000 0.000 0.000 0.000 0.000 0.000 0.000
26 1.000 0.000 0.000 0.000 0.000 0.000 0.000
27 1.000 0.000 0.000 0.000 0.000 0.000 0.000
28 1.000 0.000 0.000 0.000 0.000 0.000 0.000

# Plotly doesn't work will with both color and linetype
ggplot(U, aes(x=day, y=phome, color=trt, linetype=type)) + geom_line() +
    xl + yl + labs(color='', linetype='')

Everything has gone haywire. This is caused by the fact that both gap and day needed to be in the model to have the right main effects and interactions, but these two variables are collinear. We are requesting predictions above for combinations of day and gap that did not occur in the data, i.e., gap=1 when day> 4. Collinearities cause difficulties when geting predictions for out-of-sample data in which the collinearities disagree with those in the sample used in fitting the model.

So when the design is such that one needs to interact gap with previous state and gap is not randomly assigned or is otherwise tangled with absolute time, it is not possible to get estimates for unobserved times.

Before giving up, we’ll try a model below that does not contain gap, and that flexibly models the day main effect using a spline,

dfull <- dall[y0 == 4, .(day, y, trt, id)]
# Compute days since previous measurements
s <- dfull[, .(phome = mean(y >= 6)), by=.(trt, day)]
s[, type := 'Empirical']
empirical <- s

# It would be nice to just have day in the model and not gap (which can be
# derived from day), but the collinearities between yprev and day prevent
# this, and we must interact day with yprev - see the crosstab below

d[, table(day, yprev)]

    yprev
day    2   3   4   5   6   7
  1   32  55 224 168   0   0
  2   48  70 209 142  10   0
  3   52  73 189 113  51   0
  4   63  57 157  95 103   0
  7   69  50 132  75 148   0
  14  60  35  71  53 129 113
  28  42  12  40  37 169 147

# Treat day as numeric but with all needed variables added (less one)
# to make it be modeled as categorical except for the interaction term
# Using records truncated at death
fplain  <- lrm(y ~ gap * yprev + pol(day, 2) + trt, data=d)
# The following has a slightly greater LR chi-square:
fplain2 <- lrm(y ~ rcs(day, 4) + yprev + day %ia% yprev + trt, data=d)
# lrm(y ~ sqrt(day) * yprev + trt, data=d) also fits well
fplain

 lrm(formula = y ~ gap * yprev + pol(day, 2) + trt, data = d)
 

   1   2   3   4   5   6   7 
  50 357 297 817 526 757 489 
 

anova(fplain)

dat <- data.frame(yprev='4')
  
times <- c(1, 2, 3, 4, 7, 14, 28)
# gaps  <- c(1, diff(times))
tlev  <- levels(d$trt)

for(trt in tlev) {   # separately for treatments
    dat$trt <- trt
    dat$gap <- 1
    u <- uprob(fplain2, dat, times, 1:7, gap=TRUE)
  cat('\n', trt, '\n\n')
  print(round(u, 3))
  s <- rbind(s, data.frame(trt, day=times, p=u, phome=u[, 6] + u[, 7],
                                                 type='Model'), fill=TRUE)
}

 placebo 

       1     2     3     4     5     6     7
1  0.001 0.036 0.168 0.643 0.124 0.027 0.002
2  0.010 0.092 0.161 0.485 0.164 0.078 0.009
3  0.024 0.113 0.132 0.387 0.184 0.137 0.024
4  0.038 0.110 0.104 0.309 0.192 0.200 0.048
7  0.048 0.087 0.071 0.210 0.194 0.284 0.104
14 0.054 0.059 0.039 0.096 0.135 0.365 0.253
28 0.061 0.042 0.014 0.037 0.053 0.339 0.454

 hydroxychloroquine 

       1     2     3     4     5     6     7
1  0.001 0.037 0.169 0.642 0.123 0.026 0.002
2  0.010 0.093 0.163 0.485 0.163 0.077 0.009
3  0.025 0.114 0.133 0.387 0.182 0.135 0.024
4  0.039 0.112 0.105 0.310 0.191 0.197 0.047
7  0.050 0.089 0.072 0.212 0.194 0.282 0.102
14 0.056 0.060 0.039 0.097 0.136 0.364 0.248
28 0.063 0.043 0.015 0.037 0.054 0.341 0.448

ggplot(s, aes(x=day, y=phome, color=trt, linetype=type)) + geom_line() + xl + yl +
    labs(color='', linetype='')

for(trt in tlev) {   # separately for treatments
    dat$trt <- trt
    u <- uprob(fplain, dat, 1:28, 1:7, gap=TRUE)
  cat('\n', trt, '\n\n')
  print(round(u, 3))
  s <- rbind(s, data.frame(trt=paste(trt, 'days 1-28'), day=1:28, p=u,
                                                 phome=u[, 6] + u[, 7],
                                                 type='Model', times='Days 1-28'), fill=TRUE)
}

 placebo 

       1     2     3     4     5     6     7
1  0.001 0.033 0.155 0.643 0.136 0.030 0.002
2  0.009 0.086 0.154 0.483 0.173 0.085 0.010
3  0.022 0.108 0.131 0.386 0.185 0.142 0.026
4  0.035 0.107 0.108 0.314 0.187 0.200 0.049
5  0.046 0.094 0.088 0.254 0.181 0.255 0.082
6  0.053 0.076 0.071 0.204 0.169 0.303 0.124
7  0.058 0.056 0.055 0.160 0.152 0.341 0.177
8  0.060 0.039 0.042 0.122 0.132 0.366 0.238
9  0.062 0.025 0.029 0.091 0.111 0.377 0.305
10 0.063 0.015 0.020 0.065 0.090 0.374 0.373
11 0.064 0.008 0.012 0.045 0.071 0.361 0.439
12 0.064 0.004 0.007 0.030 0.055 0.340 0.500
13 0.064 0.002 0.004 0.019 0.043 0.315 0.552
14 0.064 0.001 0.002 0.013 0.034 0.290 0.597
15 0.064 0.000 0.001 0.009 0.027 0.266 0.633
16 0.064 0.000 0.001 0.006 0.022 0.244 0.663
17 0.064 0.000 0.000 0.005 0.019 0.224 0.688
18 0.064 0.000 0.000 0.004 0.016 0.208 0.708
19 0.064 0.000 0.000 0.003 0.015 0.194 0.724
20 0.064 0.000 0.000 0.003 0.013 0.182 0.738
21 0.064 0.000 0.000 0.003 0.012 0.172 0.749
22 0.064 0.000 0.000 0.002 0.011 0.164 0.758
23 0.064 0.000 0.000 0.002 0.011 0.158 0.765
24 0.064 0.000 0.000 0.002 0.010 0.153 0.770
25 0.064 0.000 0.000 0.002 0.010 0.150 0.774
26 0.064 0.000 0.000 0.002 0.010 0.148 0.776
27 0.064 0.000 0.000 0.002 0.010 0.147 0.777
28 0.064 0.000 0.000 0.002 0.010 0.147 0.777

 hydroxychloroquine 

       1     2     3     4     5     6     7
1  0.001 0.033 0.157 0.643 0.135 0.030 0.002
2  0.009 0.087 0.156 0.483 0.171 0.083 0.010
3  0.023 0.110 0.132 0.386 0.183 0.140 0.025
4  0.036 0.109 0.109 0.315 0.186 0.197 0.048
5  0.047 0.096 0.089 0.256 0.181 0.251 0.080
6  0.054 0.078 0.072 0.206 0.169 0.300 0.121
7  0.059 0.058 0.057 0.162 0.153 0.338 0.173
8  0.062 0.040 0.043 0.124 0.133 0.364 0.234
9  0.064 0.026 0.030 0.093 0.112 0.376 0.300
10 0.065 0.016 0.020 0.066 0.091 0.374 0.368
11 0.065 0.009 0.013 0.046 0.072 0.361 0.433
12 0.066 0.005 0.007 0.031 0.056 0.341 0.494
13 0.066 0.002 0.004 0.020 0.044 0.317 0.547
14 0.066 0.001 0.002 0.013 0.034 0.292 0.592
15 0.066 0.000 0.001 0.009 0.027 0.268 0.629
16 0.066 0.000 0.001 0.006 0.022 0.246 0.659
17 0.066 0.000 0.000 0.005 0.019 0.226 0.684
18 0.066 0.000 0.000 0.004 0.017 0.210 0.704
19 0.066 0.000 0.000 0.003 0.015 0.195 0.720
20 0.066 0.000 0.000 0.003 0.013 0.184 0.734
21 0.066 0.000 0.000 0.003 0.012 0.174 0.745
22 0.066 0.000 0.000 0.002 0.012 0.166 0.754
23 0.066 0.000 0.000 0.002 0.011 0.159 0.761
24 0.066 0.000 0.000 0.002 0.010 0.154 0.767
25 0.066 0.000 0.000 0.002 0.010 0.151 0.771
26 0.066 0.000 0.000 0.002 0.010 0.149 0.773
27 0.066 0.000 0.000 0.002 0.010 0.148 0.774
28 0.066 0.000 0.000 0.002 0.010 0.149 0.773

ggplot(s, aes(x=day, y=phome, color=trt, linetype=type)) + geom_line() + xl + yl +
    labs(color='', linetype='')

dat <- data.frame(yprev='4')
  
U <- empirical
for(trt in tlev) {   # separately for treatments
    dat$trt <- trt
    u <- uprob(fplain2, dat, 1:28, 1:7, gap=FALSE)
  U <- rbind(U, data.frame(trt, day=1:28, p=u, phome=u[, 6] + u[, 7],
                                                 type='Model days 1-28'), fill=TRUE)
}
ggplot(U, aes(x=day, y=phome, color=trt, linetype=type)) + geom_line() +
    xl + yl + labs(color='', linetype='')

Checking Validity of Method for Unequal Time Spacing Using VIOLET 2

We saw above that it is not valid to use a simple model fitted on sparse times to estimate state occupancy probabilities at unobserved times. Let’s now check that estimates are valid when observation timings are irregular but estimates are requested at actually occurring times. For this purpose we use the VIOLET 2 trial, which had a 4-level ordinal outcome scale measured on days 1-28 (records terminate at death). We’ll fit a non-covariate-adjusted (except for treatment) model using only days 1 2 3 4 7 14 28. Then we will request state occupancy probabilities at only those days, then at all days. But first start with fitting using all days and estimating all days.

Something strange happened with regard to ventilator use on day 28, so restrict analysis to first 27 days.

dcf <- readRDS('dcf.rds')
setDT(dcf)
dcf <- dcf[day < 28]

# Compute empirical proportions of being at home for those starting In Hospital, not on vent
dfull <- dcf[, .(id, day, bstatus, status, treatment, ddeath)]
emp <- dfull[bstatus == 'In Hospital',
             .(phome = mean(status == 'Home', na.rm=TRUE)),
             by=.(treatment, day)]
emp[, type := 'Empirical']
emp[, day  := unclass(day)]

dv <- dfull[day <= ddeath, ]   # don't carry death forward
s <- dv$bstatus
levels(s) <- list('Vent/ARDS'='ARDS', 'Vent/ARDS'='On Vent',
                                    'In Hospital/Facility'='In Hospital')
dv[, table(bstatus, s)]

             s
bstatus       Vent/ARDS In Hospital/Facility
  ARDS             2062                    0
  On Vent          8828                    0
  In Hospital         0                22229

dv[, bstatus := s]
setkey(dv, id, day)
dv[, pstatus := shift(status), by=id]
dv[day == 1, pstatus := bstatus]
dv[, pstatus := pstatus[drop=TRUE]]   # drop unused previous status Dead
ggplot(dv, aes(x=day)) + geom_histogram(binwidth=1) + facet_wrap(~ pstatus)

dd <- datadist(dv); options(datadist='dd')

# d2 <- dv[bstatus == 'In Hospital/Facility', ]
# Find model with best AIC
for(k in 3 : 10) {
  h <- lrm(status ~ rcs(day, k) * pstatus + treatment, data=dv, maxit=40, tol=1e-14)
  if(k == 7) f27 <- h
  cat(k, 'knots\n')
  print(AIC(h))
}

3 knots
[1] 12877.46
4 knots
[1] 12842.47
5 knots
[1] 12832.9
6 knots
[1] 12826.01
7 knots
[1] 12819.5
8 knots
[1] 12823.38
9 knots
[1] 12822.33
10 knots
[1] 12823.41

# Use 7 knots
dat <- data.frame(pstatus='In Hospital/Facility')
  
#times <- c(1, 2, 3, 4, 7, 14, 28)
tlev  <- unique(dv$treatment)
ylev  <- levels(dv$status)
times <- 1:28

R <- emp
for(trt in tlev) {   # separately for treatments
    dat$treatment <- trt
    u <- uprob(f27, dat, times, ylev, pvarname='pstatus')
  cat('\n', trt, '\n\n')
  print(round(u, 3))
  w <- data.table(treatment=trt, day=times, p=u, phome=u[, 4], type='Model')
  R <- rbind(R, w, fill=TRUE)
}

 Placebo 

    Dead Vent/ARDS In Hospital/Facility  Home
1  0.000     0.018                0.925 0.056
2  0.001     0.027                0.846 0.127
3  0.002     0.030                0.757 0.211
4  0.003     0.031                0.663 0.304
5  0.004     0.030                0.574 0.392
6  0.005     0.029                0.499 0.467
7  0.007     0.029                0.440 0.525
8  0.008     0.029                0.394 0.569
9  0.010     0.029                0.359 0.603
10 0.012     0.029                0.328 0.631
11 0.014     0.029                0.303 0.655
12 0.016     0.029                0.280 0.676
13 0.018     0.028                0.259 0.695
14 0.020     0.028                0.241 0.711
15 0.022     0.028                0.225 0.725
16 0.024     0.028                0.212 0.737
17 0.026     0.028                0.199 0.747
18 0.028     0.029                0.188 0.755
19 0.030     0.029                0.178 0.763
20 0.032     0.029                0.169 0.770
21 0.034     0.028                0.161 0.776
22 0.036     0.028                0.154 0.782
23 0.038     0.027                0.147 0.787
24 0.040     0.027                0.141 0.792
25 0.042     0.026                0.135 0.797
26 0.045     0.025                0.129 0.801
27 0.047     0.025                0.123 0.805
28 0.049     0.024                0.117 0.809

 Vitamin D 

    Dead Vent/ARDS In Hospital/Facility  Home
1  0.000     0.018                0.925 0.056
2  0.001     0.027                0.845 0.127
3  0.002     0.030                0.756 0.212
4  0.003     0.030                0.662 0.305
5  0.004     0.030                0.573 0.394
6  0.005     0.029                0.497 0.468
7  0.007     0.028                0.438 0.527
8  0.008     0.028                0.393 0.570
9  0.010     0.029                0.357 0.604
10 0.012     0.029                0.327 0.632
11 0.014     0.029                0.301 0.657
12 0.016     0.028                0.278 0.678
13 0.018     0.028                0.258 0.697
14 0.019     0.028                0.240 0.713
15 0.021     0.027                0.224 0.727
16 0.023     0.028                0.210 0.739
17 0.025     0.028                0.198 0.749
18 0.028     0.028                0.187 0.757
19 0.030     0.029                0.177 0.765
20 0.032     0.028                0.168 0.772
21 0.034     0.028                0.160 0.778
22 0.036     0.028                0.153 0.784
23 0.038     0.027                0.146 0.789
24 0.040     0.026                0.140 0.794
25 0.042     0.026                0.134 0.798
26 0.044     0.025                0.128 0.803
27 0.046     0.025                0.122 0.807
28 0.048     0.024                0.116 0.811

ggplot(R, aes(x=day, y=phome, color=treatment, linetype=type)) + geom_line() +
    xl + yl + labs(color='', linetype='')
# Fit is barely better beyond day 17 if ran lrm on data table d2

The model would probably fit better had a partial proportional odds model been used to relax the proportional odds assumption for day.

Now fit a model as if data collected was restricted to certain days. Need to recompute lags to not let a lag be based on a data that’s no longer sampled.

times <- c(1, 2, 3, 4, 7, 14, 27)
# Compute empirical proportions of being at home for those starting In Hospital, not on vent
dfull <- dcf[day %in% times, .(id, day, bstatus, status, treatment, ddeath)]
emps <- dfull[bstatus == 'In Hospital',
                            .(phome = mean(status == 'Home', na.rm=TRUE)),
                 by=.(treatment, day)]
emps[, type := 'Empirical on selected days']
emps[, day  := unclass(day)]

dv <- dfull[day <= ddeath, ]   # don't carry death forward
s  <- dv[, bstatus]
levels(s) <- list('Vent/ARDS'='ARDS', 'Vent/ARDS'='On Vent',
                                    'In Hospital/Facility'='In Hospital')
dv[, table(bstatus, s)]

             s
bstatus       Vent/ARDS In Hospital/Facility
  ARDS              611                    0
  On Vent          2408                    0
  In Hospital         0                 5925

dv[, bstatus := s]
setkey(dv, id, day)
dv[, pstatus := shift(status), by=id]
dv[day == 1, pstatus := bstatus]
dv[, pstatus := pstatus[drop=TRUE]]   # drop unused previous status Dead
dv[, gap := ifelse(day == 1, 1, day - shift(day)), by=id]
dv[, table(day, gap)]

    gap
day     1    3    7   13
  1  1352    0    0    0
  2  1336    0    0    0
  3  1321    0    0    0
  4  1303    0    0    0
  7     0 1268    0    0
  14    0    0 1212    0
  27    0    0    0 1159

ggplot(dv, aes(x=day)) + geom_histogram(binwidth=1) + facet_wrap(~ pstatus)

dd <- datadist(dv); options(datadist='dd')

# d2 <- dv[bstatus == 'In Hospital/Facility', ]
# Find transformation giving the best AIC
# Can't use a spline because of too many tied times
for(k in c(0 : 3)) {
    cat('Transformation', k, '\n')
    form <- if(k == 0) status ~ (day + gap) * pstatus + treatment
    else
        if(k == 1) status ~ (sqrt(day) + sqrt(gap)) * pstatus + treatment
    else
        if(k == 2) status ~ (pol(day, 2) + pol(gap, 2)) * pstatus + treatment
    else
        status ~ (pol(sqrt(day), 2) + pol(sqrt(day), 2)) * pstatus + treatment
    h <- lrm(form,  data=dv, maxit=50, tol=1e-14)
    if(k == 2) fsel <- h   # choose quadratic
    print(AIC(h))
}

Transformation 0 
[1] 8098.67
Transformation 1 
[1] 8008.668
Transformation 2 
[1] 7864.54
Transformation 3 
[1] 7896.565

dat <- data.frame(pstatus='In Hospital/Facility')
  
tlev  <- unique(dv$treatment)
ylev  <- levels(dv$status)

R <- rbind(R, emps, fill=TRUE)
for(trt in tlev) {   # separately for treatments
    dat$treatment <- trt
    u <- uprob(fsel, dat, times, ylev, pvarname='pstatus', gap=TRUE)
  cat('\n', trt, '\n\n')
  print(round(u, 3))
  w <- data.table(treatment=trt, day=times, p=u, phome=u[, 4], type='Model selected days')
  R <- rbind(R, w, fill=TRUE)
}

 Placebo 

    Dead Vent/ARDS In Hospital/Facility  Home
1  0.000     0.039                0.884 0.077
2  0.003     0.059                0.781 0.157
3  0.006     0.066                0.684 0.245
4  0.009     0.064                0.585 0.342
7  0.011     0.051                0.412 0.526
14 0.012     0.029                0.253 0.705
27 0.013     0.025                0.172 0.790

 Vitamin D 

    Dead Vent/ARDS In Hospital/Facility  Home
1  0.000     0.038                0.883 0.079
2  0.002     0.057                0.779 0.161
3  0.005     0.064                0.680 0.251
4  0.009     0.062                0.580 0.349
7  0.011     0.049                0.405 0.536
14 0.011     0.028                0.246 0.714
27 0.013     0.024                0.166 0.798

times <- 1:28
dat <- data.frame(pstatus='In Hospital/Facility')
for(trt in tlev) {   # separately for treatments
    dat$treatment <- trt
    u <- uprob(fsel, dat, times, ylev, pvarname='pstatus', gap=TRUE)
  cat('\n', trt, '\n\n')
  print(round(u, 3))
  w <- data.table(treatment=trt, day=times, p=u, phome=u[, 4], type='Model selected days, predict all')
  R <- rbind(R, w, fill=TRUE)
}

 Placebo 

    Dead Vent/ARDS In Hospital/Facility  Home
1  0.000     0.039                0.884 0.077
2  0.003     0.059                0.781 0.157
3  0.006     0.066                0.684 0.245
4  0.009     0.064                0.585 0.342
5  0.013     0.069                0.774 0.145
6  0.051     0.168                0.585 0.195
7  0.258     0.141                0.402 0.199
8  0.471     0.116                0.231 0.181
9  0.668     0.093                0.104 0.134
10 0.819     0.072                0.035 0.074
11 0.911     0.052                0.009 0.028
12 0.958     0.032                0.002 0.008
13 0.981     0.017                0.000 0.002
14 0.993     0.006                0.000 0.000
15 0.998     0.002                0.000 0.000
16 1.000     0.000                0.000 0.000
17 1.000     0.000                0.000 0.000
18 1.000     0.000                0.000 0.000
19 1.000     0.000                0.000 0.000
20 1.000     0.000                0.000 0.000
21 1.000     0.000                0.000 0.000
22 1.000     0.000                0.000 0.000
23 1.000     0.000                0.000 0.000
24 1.000     0.000                0.000 0.000
25 1.000     0.000                0.000 0.000
26 1.000     0.000                0.000 0.000
27 1.000     0.000                0.000 0.000
28 1.000     0.000                0.000 0.000

 Vitamin D 

    Dead Vent/ARDS In Hospital/Facility  Home
1  0.000     0.038                0.883 0.079
2  0.002     0.057                0.779 0.161
3  0.005     0.064                0.680 0.251
4  0.009     0.062                0.580 0.349
5  0.012     0.067                0.773 0.148
6  0.050     0.169                0.581 0.199
7  0.261     0.142                0.396 0.201
8  0.476     0.116                0.226 0.182
9  0.673     0.094                0.101 0.133
10 0.822     0.073                0.033 0.072
11 0.912     0.052                0.008 0.027
12 0.957     0.033                0.002 0.007
13 0.981     0.017                0.001 0.002
14 0.993     0.007                0.000 0.000
15 0.998     0.002                0.000 0.000
16 1.000     0.000                0.000 0.000
17 1.000     0.000                0.000 0.000
18 1.000     0.000                0.000 0.000
19 1.000     0.000                0.000 0.000
20 1.000     0.000                0.000 0.000
21 1.000     0.000                0.000 0.000
22 1.000     0.000                0.000 0.000
23 1.000     0.000                0.000 0.000
24 1.000     0.000                0.000 0.000
25 1.000     0.000                0.000 0.000
26 1.000     0.000                0.000 0.000
27 1.000     0.000                0.000 0.000
28 1.000     0.000                0.000 0.000

ggplot(R, aes(x=day, y=phome, color=treatment, linetype=type)) + geom_line() +
    xl + yl + labs(color='', linetype='')

This confirms that the Markov model likely handles uneven time spacings, but predictions should be limited to days appearing in the data. As before the problem is due to our attempt to make completely out-of-sample predictions on days not in the data when (1) both day and gap are important but (2) day and gap are collinear.

Better Modeling of Gap Effect in VIOLET 2

The completeness of VIOLET daily follow-up allows us to do another natural experiment where we create a full distribution of gaps by randomly sampling days within patients. For each patient randomly select a day t between 2 and L = last day of data - 6 days. Generate a random gap time r between 1 and L - 2. Remove days t, t+1, … t+r-1.

w <- dcf
w$gap <- NULL

setkey(w, id, day)
w <- w[!is.na(status), ]
# Recode baseline status like follow-up status
levels(w$bstatus) <- list('Vent/ARDS'='ARDS', 'Vent/ARDS'='On Vent',
                                            'In Hospital/Facility'='In Hospital')
w[, .(died=any(status=='Dead')), by=id][, table(died)]

died
FALSE  TRUE 
 1156   195 

u <- data.table(id=unique(w$id), key='id')
n <- nrow(u)
set.seed(1)
u[, start := sample(2 : 20, n, replace=TRUE)]
u[, width := sample(1 : 15, n, replace=TRUE)]
u[, end   := start + width - 1]
u[end < 27, table(start, width)] 

     width
start  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
   2   4  8  3  1  5  5  6  3  5  6  3  7  6  5  2
   3   6  4  4  7  6  5  4  3  6  2  6  6 13  3  2
   4   5  6  6  3  9  3  3  5  4  3  6  8  4  2  2
   5   6  4  5  4  4  6  5  5  7  4  1  3  8  3  3
   6   3  6  3  0  4  2  3  9  7  7  8  8  9  2  4
   7   1  3  5  2  3  6  7  5  6  7  5  4  7  3  4
   8   2  9  4  8  3  3  3  5  6  1  2  5  6  9  5
   9   3  3  4  5  8  3  5  4  8  5  7  3  7  4  5
   10 10  6  5  4  4  6  2  4  5  7  5  5  4  7  8
   11  4  7  4  5  6  7 10  6  6  9  3  3  4  2  1
   12  5  2  6  2  4  3  8  1  3  6  4  8  4  6  4
   13  6  3  6  4  6  4  2  4  3  8  2  7  7  2  0
   14  3  7  7  7  4  7  3  4  5  5  5  4  4  0  0
   15  8  4  8  8  7  2  7  4  4  5  4  4  0  0  0
   16  5  3  3  5  6  2  2  6  4  5  2  0  0  0  0
   17  3  9  4  6  7  1  5  4  2 10  0  0  0  0  0
   18  1  4  7  5  4  6  5  2  5  0  0  0  0  0  0
   19  3  6  4  2  0  3  4  3  0  0  0  0  0  0  0
   20  6  4  6  4  3  4  9  0  0  0  0  0  0  0  0

z <- w[u]
w <- z[end >= 27 | (day < start | day > end),]
# Recreate previous state and gap with new gaps
w[, pstatus := ifelse(day == 1, as.character(bstatus), shift(as.character(status))), by=id]
w[, gap     := ifelse(day == 1, 1L, day - shift(day)), by=id]
# Don't carry death forward 
w[,  nd := cumsum(status == 'Dead'), by=id]
w <- w[nd < 2, ]
w[, sum(nd > 0)]    # original number of deaths

[1] 195

w[day > 1, table(gap)]

gap
    1     2     3     4     5     6     7     8     9    10    11    12    13 
22788    76    92    84    77    89    73    85    73    78    84    61    71 
   14    15    16 
   78    45    36 

with(w, ggfreqScatter(day, gap))
dd <- datadist(w); options(datadist='dd')
f0 <- lrm(status ~ pstatus + rcs(day, 7) + treatment + rcs(age, 4) + sofa + lips,
         data=w, tol=1e-13)
f1 <- lrm(status ~ rcs(gap, 4) + pstatus + rcs(day, 4) + treatment + rcs(age, 4) + sofa + lips,
         data=w, tol=1e-13)
f2 <- lrm(status ~ gap * pstatus + rcs(day, 4) + treatment + rcs(age, 4) + sofa + lips,
         data=w, tol=1e-13)
f3 <- lrm(status ~ lsp(gap, 2) * pstatus + rcs(day, 4) + treatment + rcs(age, 4) + sofa + lips,
         data=w, tol=1e-13)
f4 <- lrm(status ~ lsp(gap, c(2,4)) * pstatus + rcs(day, 4) + treatment + rcs(age, 4) + sofa + lips,
         data=w, tol=1e-13)
f5 <- lrm(status ~ rcs(gap, 4) * pstatus + rcs(day, 4) + treatment + rcs(age, 4) + sofa + lips,
         data=w, tol=1e-13)
c(f0=AIC(f0), f1=AIC(f1), f2=AIC(f2), f3=AIC(f3), f4=AIC(f4), f5=AIC(f5))

      f0       f1       f2       f3       f4       f5 
12012.40 11652.23 11625.38 11597.44 11601.45 11600.24 

ggplot(Predict(f1, gap, conf.int=FALSE))
ggplot(Predict(f5, gap, pstatus, conf.int=FALSE))

anova(f5)

Before drawing conclusions let’s handle non-PO in time.

require(VGAM)

fpol <- vglm(ordered(status) ~ pstatus + day + treatment + rcs(age, 4) + sofa + lips,
          cumulative(reverse=TRUE, parallel=TRUE), data=w)
fpo3 <- vglm(ordered(status) ~ pstatus + rcs(day, 3) + treatment + rcs(age, 4) + sofa + lips,
          cumulative(reverse=TRUE, parallel=TRUE), data=w)
fpo5 <- vglm(ordered(status) ~ pstatus + rcs(day, 5) + treatment + rcs(age, 4) + sofa + lips,
          cumulative(reverse=TRUE, parallel=TRUE), data=w)
f0l <- vglm(ordered(status) ~ pstatus + day + treatment + rcs(age, 4) + sofa + lips,
          cumulative(reverse=TRUE, parallel=FALSE ~ day), data=w)
f03 <- vglm(ordered(status) ~ pstatus + rcs(day, 3) + treatment + rcs(age, 4) + sofa + lips,
          cumulative(reverse=TRUE, parallel=FALSE ~ rcs(day, 3)), data=w) 
f05 <- vglm(ordered(status) ~ pstatus + rcs(day, 5) + treatment + rcs(age, 4) + sofa + lips,
          cumulative(reverse=TRUE, parallel=FALSE ~ rcs(day, 5)), data=w)
f1 <- vglm(ordered(status) ~ pstatus + rcs(gap, 4) + rcs(day, 3) + treatment + rcs(age, 4) + sofa + lips,
          cumulative(reverse=TRUE, parallel=FALSE ~ rcs(day, 3)), data=w)
f2 <- vglm(ordered(status) ~ gap * pstatus + rcs(day, 3) + treatment + rcs(age, 4) + sofa + lips,
          cumulative(reverse=TRUE, parallel=FALSE ~ rcs(day, 3)), data=w)
f3 <- vglm(ordered(status) ~ rcs(gap, 3) * pstatus + rcs(day, 3) + treatment + rcs(age, 4) + sofa + lips,
          cumulative(reverse=TRUE, parallel=FALSE ~ rcs(day, 3)), data=w)
f4 <- vglm(ordered(status) ~ rcs(gap, 4) * pstatus + rcs(day, 3) + treatment + rcs(age, 4) + sofa + lips,
          cumulative(reverse=TRUE, parallel=FALSE ~ rcs(day, 3)), data=w)

c(fpol=AIC(fpol), fpo3=AIC(fpo3), fpo5=AIC(fpo5), f0l=AIC(f0l), f03=AIC(f03), f05=AIC(f05), f1=AIC(f1), f2=AIC(f2), f3=AIC(f3), f4=AIC(f4))

    fpol     fpo3     fpo5      f0l      f03      f05       f1       f2 
12044.31 12023.26 12023.96 12039.35 11988.62 11992.68 11639.42 11612.38 
      f3       f4 
11590.84 11592.26 

Comparing AICs indicates the following:

• fpol vs. fpo3: effect of day is nonlinear in a PO model
• fpo3 vs. fpo5: 3 knots are adequate for day in PO model
• fpo5 vs. f05: moderate non-PO for day
• f05 vs. f0l: non-PO effect of day or main effect of day or both are nonlinear
• f05 vs. f03: 3 knots works better for day for combined main and non-PO effects
• f1 vs. f03: main effect is needed for gap
• f2 vs. f1: previous state interacts weekly with gap
• f2, f3,f4: previous state interaction with gap is better represented by 3 knots

Odds Ratios for ORCHID

Now get treatment odds ratios for state occupancy probabilities, as a function of time, as opposed to odds of transitions over consecutive time periods.

odds <- function(x) x / (1. - x)
setDT(Rorchid, key=c('trt', 'day'))
b    <- Rorchid[trt == 'hydroxychloroquine', phome]
a    <- Rorchid[trt == 'placebo',            phome]
days <- Rorchid[trt == 'placebo',            day]
or   <- odds(b) / odds(a)
plot(days, or, type='b', xlab='Day', ylab='hydroxy : placebo OR', ylim=c(1, 1.05))

abline(h=exp(coef(f)['trt=hydroxychloroquine']), col=gray(0.85))
text(27, 1.007,
         'Horizontal line is conditional\ntransition OR from unified PO model', adj=1)

Bayesian Markov Model for ORCHID

We parallel the regular proportional odds model with a Bayesian Markov proportional hazards model.

stanSet()

bmarkov <- blrm(y ~ gap * yprev + pol(day, 2) + trt + rcs(age, 3) +
                sofa_nogcs,
                data=d, refresh=100, file='bmarkov.rds')
htmlVerbatim(stanDx(bmarkov))

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>=2                    2662 1.001
 y>=3                    2219 1.001
 y>=4                    2691 1.000
 y>=5                    3632 1.001
 y>=6                    3459 1.001
 y>=7                    2280 1.001
 gap                     2625 1.000
 yprev=3                 3087 1.000
 yprev=4                 3970 1.000
 yprev=5                 3637 1.001
 yprev=6                 2769 1.001
 yprev=7                 2578 1.000
 day                     5192 1.000
 day^2                   5053 1.000
 trt=hydroxychloroquine  5442 1.000
 age                     5505 1.000
 age'                    5513 1.000
 sofa_nogcs              4865 1.000
 gap * yprev=3           4791 1.000
 gap * yprev=4           4223 1.000
 gap * yprev=5           5252 0.999
 gap * yprev=6           5020 0.999
 gap * yprev=7           5482 0.999
 

stanDxplot(bmarkov)
bmarkov

Bayesian Proportional Odds Ordinal Logistic Model

Dirichlet Priors With Concentration Parameter 0.308 for Intercepts

 blrm(formula = y ~ gap * yprev + pol(day, 2) + trt + rcs(age, 
     3) + sofa_nogcs, data = d, refresh = 100, file = "bmarkov.rds")
 

Frequencies of Responses

   1   2   3   4   5   6   7 
  50 357 297 817 526 757 489 
 

We need to be able to compute posterior distributions on state occupancy probabilities to take into account all sources of uncertainty. The follow function computes all needed derived quantities, separately for each posterior draw.

# Compute unconditional probabilities
# Output is an array of dimension number of draws by s rows corresponding to
# t=times(1), times(2), ...
# times(s), and final dimension is the number of Y categories
# times(1) must be 1
# 
# fit is the Bayesian model fit from blrm which must contain the variable named day and the
# previous state, which must be named by the value of pvarname
# data is a data frame or data table that contains covariate settings
# other than day and what's in pvarname.  data has one row.
# times is the vector of days for which to evaluate state occupancy probabilities
# ylevels are the ordered levels of the outcome variable
# Set gap=TRUE to compute time gaps when predictor gap is in the model

uprobB <- function(fit, data, times, ylevels, pvarname='yprev', gap=FALSE) {
    if(! length(fit$draws)) stop('fit does not have posterior draws')
    nd <- nrow(fit$draws)
  k  <- length(ylevels)
  s  <- length(times)
  P  <- array(NA, c(nd, s, k),
                        dimnames=list(paste('draw', 1 : nd), as.character(times), 
                                                    as.character(ylevels)))
  # Never uncondition on initial state
  data$day <- times[1]
  if(gap) data$gap <- times[1]
  p <- predict(fit, data, type='fitted.ind', posterior.summary='all')
  P[, 1, ] <- p
  # cp: matrix of conditional probabilities of Y conditioning on previous time Y
  # Columns = k conditional probabilities conditional on a single previous state
  # Rows    = all possible previous states
  # This is for a single posterior draw
  cp <- matrix(NA, nrow=k, ncol=k, 
               dimnames=list(paste('t-1', ylevels), paste('t', ylevels)))
  data <- as.list(data)
  data[[pvarname]] <- as.character(ylevels[-1])
  # don't request estimates after absorbing state
  edata <- expand.grid(data)
  for(it in 2 : s) {
    edata$day <- times[it]
    if(gap) edata$gap <- times[it] - times[it - 1]
    # y=first level = absorbing state
    pp <- predict(fit, edata, type='fitted.ind', posterior.summary='all')
    for(idraw in 1 : nd) {
        cp <- rbind(c(1., rep(0., k - 1)), pp[idraw, ,])
      # Compute unconditional probabilities of being in all possible states
      # at current time t
      P[idraw, it, ] <- t(cp) %*% P[idraw, it - 1, ]
    }
  }
  P
}

times <- c(1,2,3,4,7,14,28)
ylev <- 1 : 7
tlev <- levels(d$trt)
dat <- data.frame(yprev='4', age=repcov$age, sofa_nogcs=repcov$sofa_nogcs)
R   <- list()
for(trt in tlev) {   # separately for treatments
    dat$trt <- trt
    u <- uprobB(bmarkov, dat, times, ylev, pvarname='yprev', gap=TRUE)
    R[[trt]] <- u
}

# Compute summary measure per treatment, separately for each time and posterior draw
# P(home) = P(Y=6 or Y=7)

s <- function(u) u[, 6] + u[, 7]
# w <- apply(R$placebo, 1, s)    # to test
r <- lapply(R, function(x) apply(x, 1, s))

For each treatment and time compute the posterior mean P(home) and 0.95 highest posterior density intervals.

g <- function(x) {
    w <- c(mean(x), HPDint(x))
    names(w) <- c('mean', 'lower', 'upper')
    w
}

rs <- lapply(r, function(x) as.data.frame(t(apply(x, 1, g))))
for(trt in tlev) rs[[trt]]$trt <- trt
rs <- do.call(rbind, rs)
# Fetch day from end of row names (follows .)
rs$day <- as.numeric(sub('^.*\\.(.*?)', '\\2', rownames(rs), ))
xl <- xlab('Day'); yl <- ylab('P(being at home)')
rs$txt <- with(rs, paste0(trt, '<br>Day:', day, '<br>', round(mean, 3)))
gg <- ggplot(rs, aes(x=day, y=mean, color=trt, label=txt)) + geom_line() + 
    geom_ribbon(aes(ymin = lower, ymax = upper), alpha=0.2, linetype=0) +
    xl + yl + labs(color='')
ggp(gg)

Instead of individual per-treatment estimates summarize the difference and odds ratio for state occupancy probabilities.

Dif <- r$hydroxychloroquine - r$placebo
dif <- as.data.frame(t(apply(Dif, 1, g)))
dif$day <- as.numeric(sub('^.*\\.(.*?)', '\\2', rownames(dif), ))
dif$txt <- with(dif, paste0('Day:', day, '<br>', round(mean, 3)))
gg <- ggplot(dif, aes(x=day, y=mean, label=txt)) + geom_line() +
    geom_ribbon(aes(ymin=lower, ymax=upper), alpha=0.2, linetype=0) +
    xl + ylab('Hydroxy - Placebo Difference in P(home)')
ggp(gg)

or <- odds(r$hydroxychloroquine) / odds(r$placebo)
or <- as.data.frame(t(apply(or, 1, g)))
or$day <- as.numeric(sub('^.*\\.(.*?)', '\\2', rownames(or), ))
or$txt <- with(or, paste0('Day:', day, '<br>', round(mean, 3)))
gg <- ggplot(or, aes(x=day, y=mean, label=txt)) + geom_line() +
    geom_ribbon(aes(ymin=lower, ymax=upper), alpha=0.2, linetype=0) +
    xl + ylab('Hydroxy - Placebo OR for Being at Home')
ggp(gg)

Evidence for Efficacy Over Time for One Covariate Combination

Now compute the Bayesian evidence measure for a positive difference in P(home). Note that in the two scatterplots relating posterior draws for differences in occupancy probabilities (P(home)), plotting one day vs. another, no points are on the upper left or lower right quadrants, so the evidence for a positive treatment effect does not vary over the days of follow-up.

plot(Dif[1, ], Dif[5, ], xlab='Posterior Draws for Day 1', ylab='Day 7', cex=.3)
abline(h=0, v=0, col=gray(0.8))
plot(Dif[3, ], Dif[6, ], xlab='Posterior Draws for Day 3', ylab='Day 14', cex=.3)
abline(h=0, v=0, col=gray(0.8))
pp <- apply(Dif, 1, function(x) mean(x > 0))
w <- data.frame(day = as.integer(names(pp)), pp)
w$txt <- with(w, paste0('Day:', day, '<br>', round(pp, 3)))
gg <- ggplot(w, aes(x=day, y=pp, label=txt)) + geom_line() +
    xl + ylab('P(increase in P(home) with hydroxychloroquine)')
ggp(gg)

The evidence does vary when computing the posterior probability that treatment has a benefit exceeding an absolute risk difference of 0.02. This is primarily due to the near impossibility of being discharged on the first or second day.

pp <- apply(Dif, 1, function(x) mean(x > 0.02))
w <- data.frame(day = as.integer(names(pp)), pp)
w$txt <- with(w, paste0('Day:', day, '<br>', round(pp, 3)))
gg <- ggplot(w, aes(x=day, y=pp, label=txt)) + geom_line() +
    xl + ylab('P(increase in P(home) with hydroxychloroquine > 0.02)')
ggp(gg)

Posterior Distributions of Median and Restricted Mean Time to Home

Even though the data collection times in ORCHID were not granular enough to compute median time until the patient returned home, we compute this anyway, and compute per-treatment posterior means and discrete posterior probability distributions for these medians. For a given posterior sample, the median is the first day for which P(home) >= 0.5 for a treatment. We arbitrary use 30d as the median when the probability never exceeds 0.5. Recall that all these calculations are conditional on the patient being in state 4 (hospitalized on supplemental oxygen) pre-randomization.

days <- c(1,2,3,4,7,14,28)
g <- function(p) if(all(p < 0.5)) 30 else min(days[p >= 0.5])
rmed <- lapply(r, function(x) data.frame(med = apply(x, 2, g)))
for(trt in tlev) rmed[[trt]]$trt <- trt
rmed <- do.call(rbind, rmed)
with(rmed, tapply(med, trt, mean))   # posterior mean for median

hydroxychloroquine            placebo 
                14                 14 

with(rmed, table(med, trt))          # discrete distribution for median

    trt
med  hydroxychloroquine placebo
  14               4000    4000

The sharpness of the posterior distribution of the median time to home is just a reflection of what an easy question we are asking: what is the median of time to home after essentially rounding it to the nearest 7d.

Follow-up was not long enough to compute the mean time until returning home, but we can do inference for the mean restricted time until home (restricted to 28d) assuming that home is an absorbing state (which is approximately but not exactly true). Let \(t_{i}, i=1, 2, \dots, 7\) represent the assessment times (days 1, 2, 3, 4, 7, 14, 28) for time until home \(T\). Under the \(T \leq 28\) restriction, the mean restricted time to home is \(\sum_{i=1}^{7} t_{i} \Pr(T = t_{i}) = \sum_{i=1}^{7} t_{i} [\Pr(T \leq t_{i}) - [i > 1] \Pr(T \leq t_{i-1})]\).

g <- function(p) sum(days * (p - c(0., p[-length(p)])))
rmean <- lapply(r, function(x) data.frame(mean = apply(x, 2, g)))
for(trt in tlev) rmean[[trt]]$trt <- trt
rmean <- do.call(rbind, rmean)
with(rmean, tapply(mean, trt, mean))   # posterior means for mean restricted time to home

hydroxychloroquine            placebo 
          10.55112           10.57485 

ggplot(rmean, aes(x=mean, color=trt)) + geom_density() +
    labs(color='') + xlab('Restricted Mean Time to Home (days)') + ylab('')

Compute the posterior probability that the mean restricted time to home is shorter for hydroxy.

a <- rmean[grep('placebo', rownames(rmean)), 'mean']
b <- rmean[grep('hydrox',  rownames(rmean)), 'mean']
plot(density(b - a), main='', xlab='Difference in Mean Restricted Time to Home', ylab='')
abline(v=0, col=gray(0.85))
mean(b < a)

[1] 0.56825

Markov Model with Random Effects

Let’s see what difference it makes to add patient-level random effects to the Markov state transition model.

bmarkovre <- blrm(y ~ gap * yprev + pol(day, 2) + trt + rcs(age, 3) +
                sofa_nogcs + cluster(id),
                data=d, refresh=100, file='bmarkovre.rds')
htmlVerbatim(stanDx(bmarkovre))

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>=2                    2058 1.000
 y>=3                    1266 1.002
 y>=4                    1485 1.001
 y>=5                    2911 1.000
 y>=6                    2458 1.000
 y>=7                    1585 1.001
 gap                     1950 1.001
 yprev=3                 1880 1.000
 yprev=4                 3299 1.000
 yprev=5                 2750 1.001
 yprev=6                 1627 1.002
 yprev=7                 2037 1.000
 day                     5880 0.999
 day^2                   5076 1.001
 trt=hydroxychloroquine  5040 0.999
 age                     4770 1.001
 age'                    5632 0.999
 sofa_nogcs              5487 0.999
 gap * yprev=3           4795 1.000
 gap * yprev=4           5105 0.999
 gap * yprev=5           4534 0.999
 gap * yprev=6           5255 1.000
 gap * yprev=7           5623 0.999
 sigmag                  1384 1.000
 

stanDxplot(bmarkovre)
bmarkovre

Bayesian Proportional Odds Ordinal Logistic Model

Dirichlet Priors With Concentration Parameter 0.308 for Intercepts

 blrm(formula = y ~ gap * yprev + pol(day, 2) + trt + rcs(age, 
     3) + sofa_nogcs + cluster(id), data = d, refresh = 100, file = "bmarkovre.rds")
 

Frequencies of Responses

   1   2   3   4   5   6   7 
  50 357 297 817 526 757 489 
 

The very low standard deviation of the random effects (0.053) lends credence to the idea that the state transitions are conditionally independent and that random effects are not needed. Also, the posterior mean and median regression parameter values and standard deviations of their posterior distributions changed very little.

Serial Transition Model With Non-Proportional Odds

Return to the Markov model without random effects. Let’s use a partial proportional odds model in which we add departures from the proportional odds assumption with regard to two variables: previous state and time (but not their interaction). We are using the constrained partial proportional odds model in which the constraint imposes linearity in the outcome category (but not for the effect of the previous outcome level).

bmarkovppo <- blrm(y ~ gap * yprev + pol(day, 2) + trt + rcs(age, 3) +
                sofa_nogcs, ~ day + yprev, cppo=function(y) y,
                data=d, refresh=100, file='bmarkovppo.rds')
htmlVerbatim(stanDx(bmarkovppo))

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>=2                     524 1.009
 y>=3                     480 1.008
 y>=4                     441 1.009
 y>=5                     589 1.005
 y>=6                     795 1.001
 y>=7                     812 1.003
 gap                      487 1.010
 yprev=3                  978 1.004
 yprev=4                  595 1.006
 yprev=5                  761 1.004
 yprev=6                  902 1.005
 yprev=7                  631 1.009
 day                     2249 1.000
 day^2                   1787 1.001
 trt=hydroxychloroquine  2581 1.000
 age                     2304 1.002
 age'                    2419 1.000
 sofa_nogcs              1996 1.001
 gap * yprev=3           1883 1.003
 gap * yprev=4           2312 0.999
 gap * yprev=5           2078 1.000
 gap * yprev=6           1248 1.001
 gap * yprev=7           1702 1.001
 day x f(y)               485 1.008
 yprev=3 x f(y)           886 1.002
 yprev=4 x f(y)           519 1.008
 yprev=5 x f(y)          1056 1.003
 yprev=6 x f(y)           852 1.005
 yprev=7 x f(y)           549 1.008
 

stanDxplot(bmarkovppo)
bmarkovppo

Bayesian Constrained Partial Proportional Odds Ordinal Logistic Model

Dirichlet Priors With Concentration Parameter 0.308 for Intercepts

 blrm(formula = y ~ gap * yprev + pol(day, 2) + trt + rcs(age, 
     3) + sofa_nogcs, ppo = ~day + yprev, cppo = function(y) y, 
     data = d, refresh = 100, file = "bmarkovppo.rds")
 

Frequencies of Responses

   1   2   3   4   5   6   7 
  50 357 297 817 526 757 489 
 

# htmlVerbatim(compareBmods(bmarkovppo, bmarkov))
# Note: bmarkovppo fit with loo=TRUE was not saved because it was 171MB for some reason
# Output:
# Method: stacking
# ------
#        weight
#  model1 0.957 
#  model2 0.043 
#  
#  Compute posterior means divided by SE
round(coef(bmarkovppo) / sqrt(diag(vcov(bmarkovppo))), 2)

              y>=2                   y>=3                   y>=4 
              8.57                  -4.34                  -9.74 
              y>=5                   y>=6                   y>=7 
            -16.33                 -17.26                 -18.16 
               gap                yprev=3                yprev=4 
              0.95                   8.07                  18.18 
           yprev=5                yprev=6                yprev=7 
             25.63                  25.42                  11.92 
               day                  day^2 trt=hydroxychloroquine 
              8.53                 -10.99                   0.20 
               age                   age'             sofa_nogcs 
             -2.07                  -0.34                  -4.84 
     gap * yprev=3          gap * yprev=4          gap * yprev=5 
             -1.64                  -4.98                  -6.07 
     gap * yprev=6          gap * yprev=7             day x f(y) 
            -10.44                  -5.24                  17.72 
    yprev=3 x f(y)         yprev=4 x f(y)         yprev=5 x f(y) 
             -0.32                   1.02                  -3.09 
    yprev=6 x f(y)         yprev=7 x f(y) 
             -6.69                  -2.70 

There were some minor sampling/convergence problems especially with respect to the parameters for departures from proportional odds.

The model comparison results listed in the code above puts a great likelihood on the partial proportional odds model fitting better than the proportional odds model. The result of compareBmods, not saved in the code because the loo=TRUE option makes the fit object much bigger, is here (but from a slightly different fit):

Method: stacking
------
       weight
 model1 0.957 
 model2 0.043 

Based on posterior mean parameter values and standard deviations of posterior distributions (standard errors), the most important departures from PO seem to be time-related and the effect of previously being in state 6.

Final Model With Non-PO With Respect to Time

Our final model will have as non-proportional odds only the effect of time, but we will now allow the time effect to be nonlinear (quadratic). As before the cppo function assumes that the PO departure is proportional to the numeric 1-7 outcome codes. A fully unconstrained partial proportional odds model would have too many parameters.

bmarkovppof <- blrm(y ~ gap * yprev + pol(day, 2) + trt + rcs(age, 3) +
                sofa_nogcs, ~ pol(day, 2), cppo=function(y) y,
                data=d, refresh=100, file='bmarkovppof.rds')
htmlVerbatim(stanDx(bmarkovppof))

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>=2                    3155 1.001
 y>=3                    2772 1.001
 y>=4                    3497 1.000
 y>=5                    4565 1.000
 y>=6                    3112 1.000
 y>=7                    2744 1.000
 gap                     3701 1.000
 yprev=3                 3236 1.001
 yprev=4                 4312 1.000
 yprev=5                 3439 1.001
 yprev=6                 2914 1.001
 yprev=7                 3443 1.000
 day                     5998 1.000
 day^2                   4539 1.000
 trt=hydroxychloroquine  7707 0.999
 age                     6236 1.000
 age'                    6518 1.000
 sofa_nogcs              5643 1.000
 gap * yprev=3           5416 0.999
 gap * yprev=4           7418 1.000
 gap * yprev=5           6227 0.999
 gap * yprev=6           4029 1.000
 gap * yprev=7           4523 1.000
 day x f(y)              2652 1.000
 day^2 x f(y)            4295 1.001
 

stanDxplot(bmarkovppof)
bmarkovppof

Bayesian Constrained Partial Proportional Odds Ordinal Logistic Model

Dirichlet Priors With Concentration Parameter 0.308 for Intercepts

 blrm(formula = y ~ gap * yprev + pol(day, 2) + trt + rcs(age, 
     3) + sofa_nogcs, ppo = ~pol(day, 2), cppo = function(y) y, 
     data = d, refresh = 100, file = "bmarkovppof.rds")
 

Frequencies of Responses

   1   2   3   4   5   6   7 
  50 357 297 817 526 757 489 
 

There is strong evidence that the departure from PO with respect to day of follow-up is nonlinear.

Final Check

Let’s fit a Bayesian transition model without covariates (other than treatment and initial state) to once again check against empirical proportions of being at home.

bmarkovppon <- blrm(y ~ gap * yprev + pol(day, 2) + trt, ~ pol(day, 2),
                                        cppo=function(y) y,
                      data=d, refresh=100, file='bmarkovppon.rds')
htmlVerbatim(stanDx(bmarkovppon))

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>=2                    3039 1.000
 y>=3                    2804 1.001
 y>=4                    3198 1.001
 y>=5                    4132 1.000
 y>=6                    3108 1.000
 y>=7                    2452 1.000
 gap                     3223 1.002
 yprev=3                 2929 1.000
 yprev=4                 3561 1.000
 yprev=5                 3773 1.000
 yprev=6                 2980 0.999
 yprev=7                 3157 1.000
 day                     5230 1.000
 day^2                   4511 0.999
 trt=hydroxychloroquine  5522 1.000
 gap * yprev=3           4918 0.999
 gap * yprev=4           5557 1.000
 gap * yprev=5           4704 1.000
 gap * yprev=6           3538 1.000
 gap * yprev=7           3973 1.001
 day x f(y)              2504 1.000
 day^2 x f(y)            4184 0.999
 

bmarkovppon

Bayesian Constrained Partial Proportional Odds Ordinal Logistic Model

Dirichlet Priors With Concentration Parameter 0.308 for Intercepts

 blrm(formula = y ~ gap * yprev + pol(day, 2) + trt, ppo = ~pol(day, 
     2), cppo = function(y) y, data = d, refresh = 100, file = "bmarkovppon.rds")
 

Frequencies of Responses

   1   2   3   4   5   6   7 
  50 357 297 817 526 757 489