Introduction

The NIH NIAID ACTT-1 study was a double-blind randomized clinical trial of IV remdesivir vs. placebo for treatment of hospitalized COVID-19 patients with evidence of lower respiratory tract infection. The main publication is here. The primary endpoint was time to recovery, defined by either discharge from the hospital or hospitalization for infection control rather than for treatment. 1062 patients were randomized (541 to remdesivir and 521 to placebo). The rate ratio for recovery was 1.29 with 0.95 CI of [1.12, 1.49], P < 0.001 by a log-rank test. In a secondary analysis using the proportional odds model with an 8-category ordinal scale assessed on day 15 the odds ratio was 1.5 with 0.95 CI of [1.2, 1.9]. The 15-day cumulative incidence of mortality was 0.067 with remdesivir and 0.119 with placebo (hazard ratio 0.73, [0.52, 1.03]).

Eleven patients with missing baseline ordinal score and all subsequent ordinal outcome assessments also missing were excluded. Analyses that followed used 1051 patients (533 on remdesivir and 518 on placebo).

Patients dishcharged to home were only assessed on certain days from randomization. Days after the last day in hospital (outcome status y > 2) up until death or end of follow-up for which no formal outcome assessment was done have records with y = 1.5, indicating at home but need for O$_2$ or assistance status is unknown. When doing a Bayesian analysis with the rmsb package blrm function these observations are considered to be interval censored in [1,2]. When a previous state is interval censored, the previous state is taken to be a new state, 1.5.

After much experimentation with various longitudinal models (and also see here and here), a first-order Markov state transition model for ordinal outcome states provides the best fit to the correlation structure found in the data. The Markov model uses the previous ordinal state as a (categorical) predictor of the current state. It allows the correlation between states on two successive days to be arbitrarily high within patient, and it elegantly handles absorbing states such as death. Excellent handling of the correlation structure is done without the complexities of random effects. Go here for a list of papers about Markov longitudinal modeling of binary and ordinal outcomes.

Markov models in their usual form assume that the outcome variable is also assessed at baseline (time zero). Such is the case for ACTT-1. The fact that some of the outcome levels (e.g., death) are not used at baseline is of no consequence.

The only obstacle to using transition models is that the primary model parameters (here, in a proportional odds model) relate to the probability of being in the current state or higher given the previous state. For day 1 post randomization, this is no problem. But we need to restate the model to obtain time-unconditional probabilities of being in state y or worse (“state occupancy probabilities”). This is done through recursive (through time) matrix multiplications to “uncondition” the transition probabilities to get probabilities of being in the current state regardless of previous states. This are simple to interpret quantities that are functions of all model parameters. Because of that it is difficult to get confidence intervals and statistical tests using frequentist methods. In examples below it is seen that the Bayesian approach provides exact inference, taking all uncertainties into account. Graphs will show the probability of the patient being at home as a function of time and treatment. Home is almost, but not quite, an absorbing state. So it is possible that such graphs will increase up to a point then slightly decrease (that didn’t happen in this study). If we consider home to be an absorbing state, then the graph is a cumulative incidence graph.

Various ways of restating the model results are also given below. One of these is the mean restricted time to the patient going home.

Note that the transition model needs to take absolute time (days from randomization) into account, as well as previous state. That is because the mix of outcomes changes over time in a way that simple transition odds ratios can’t capture.

To be added to this report is an example of using the fitted model to simulate longitudinal ordinal data.

require(rmsb)

require(data.table)

require(plotly)

knitrSet(lang='markdown', w=7, h=7, dev='png', fig.path='png/a-')
options(prType='html')
ggp <- ggplotlyr   # ggplotlyr is in Hmisc
# Gets around a bug in tooltip construction with ggplotly
    
d <- readRDS('actt1.rds')
ordlabels <- c('not hosp, no limits', 'not hosp', 'not hosp, limits or O2',
  'hosp, no O2, no care', 'hosp, no O2, care',
  'hosp, O2', 'hosp, noninvasive vent', 'hosp, invasive vent', 'death')

# Ccreate a data table with information that doesn't change over days
u <- d[seq == 1, .(id, tx, ddeath, dlasthosp, dlasthospf, ttrecov0)]
setkey(u, 'id')

# Also create a data table that carries death forward to day 28.

w <- d[, .(dlast=day, last8 = day == max(day) & y == 8,
           tx=tx, y0=y0, dlasthosp, dlasthospf, ttrecov0),
                     by=id]
w <- w[last8 == TRUE, .(day = (dlast + 1) : 28,
                        y   = rep(8, 28 - dlast), tx=tx, y0=y0,
                                                dlasthosp=dlasthosp, dlasthospf=dlasthospf,
                                                ttrecov0=ttrecov0),
             by=id]
label(w$day) <- label(d$day)
label(w$y)   <- label(d$y)
w <- rbind(d[, .(id, day, y, tx, y0, dlasthosp, dlasthospf, ttrecov0)], w)
setkey(w, id, day)
dcf <- w

Descriptive Statistics

html(describe(d))

d

23 Variables 26167 Observations

id

n	missing	distinct
26167	0	1051

lowest :	COV.00701	COV.00702	COV.00703	COV.00704	COV.00705
highest:	COV.01802	COV.01803	COV.01805	COV.01808	COV.01809

tx: Assigned Treatment Format:$

n	missing	distinct
26167	0	2

 Value         Placebo Remdesivir
 Frequency       12743      13424
 Proportion      0.487      0.513

ddeath: Day of First Record With Status = Death

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
1674	24493	27	0.997	16.93	8.534	5	6	11	17	23	27	29

lowest : 2 3 4 5 6 , highest: 25 26 27 29 31

dhome: First day with status at home (NA if never)

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
19070	7097	32	0.974	17.25	5.53	12	13	14	15	21	27	28

lowest : 2 4 5 6 7 , highest: 30 31 32 33 40

dlasthosp: Last day in hospital with outcome assessment

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
26167	0	33	0.997	14.29	10.73	2	3	6	12	24	28	29

lowest : 0 1 2 3 4 , highest: 28 29 30 31 32

dlasthospf: Last day in first hospital admission with outcome assessment

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
26167	0	33	0.997	13.86	10.56	2	3	5	11	22	28	29

lowest : 0 1 2 3 4 , highest: 28 29 30 31 32

d29dthe0: Time to Death days

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
1646	24521	26	0.997	15.55	8.174	4	5	10	15	22	25	26

lowest : 1 2 3 4 5 , highest: 23 24 25 26 28

d29dthe1: Time to Death Censor Time Point days

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
24521	1646	23	0.504	27.45	0.9683	25	26	28	28	28	28	28

lowest : 0 1 2 3 4 , highest: 24 25 26 27 28

actarm: Treatment

n	missing	distinct
26167	0	4

 Value              Not Treated             Placebo          Remdesivir
 Frequency                    4               12017               13421
 Proportion               0.000               0.459               0.513
                               
 Value      Unplanned Treatment
 Frequency                  725
 Proportion               0.028

y0: Baseline Status

n	missing	distinct	Info	Mean	Gmd
26167	0	4	0.898	5.56	1.12

 Value          4     5     6     7
 Frequency   3651 11094  4531  6891
 Proportion 0.140 0.424 0.173 0.263

ttrecov0: Time to Recovery days

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
20486	5681	33	0.997	10.16	7.918	2	3	4	8	14	21	25

lowest : 0 1 2 3 4 , highest: 28 29 30 31 32

ttrecov1: Time to Recovery Censor Time Point days

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
5681	20486	18	0.859	27.1	1.928	25	27	27	27	28	28	30

lowest : 0 1 2 3 4 , highest: 25 26 27 28 30

 Value          0     1     2     3     4     5     6     7    10    13    14    15
 Frequency      7     6    12    10     6     6    15     8    12    14    29    17
 Proportion 0.001 0.001 0.002 0.002 0.001 0.001 0.003 0.001 0.002 0.002 0.005 0.003
                                               
 Value         24    25    26    27    28    30
 Frequency     50   156   135  2464  2203   531
 Proportion 0.009 0.027 0.024 0.434 0.388 0.093

y15: Ordinal Scale Score Day 15

n	missing	distinct	Info	Mean	Gmd
26167	0	8	0.955	3.597	2.693

lowest : 1 2 3 4 5 , highest: 4 5 6 7 8

 Value          1     2     3     4     5     6     7     8
 Frequency   7479  5934   586  1862  3033  1231  5322   720
 Proportion 0.286 0.227 0.022 0.071 0.116 0.047 0.203 0.028

age: Age

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
26167	0	68	1	57.98	16.62	32	38	48	59	68	77	82

lowest : 21 22 23 25 26 , highest: 86 87 88 89 90

sex: Sex

n	missing	distinct
26167	0	2

 Value          F     M
 Frequency   9309 16858
 Proportion 0.356 0.644

dursx: Duration of Symptoms days

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
26083	84	35	0.995	9.865	5.669	3	4	6	9	12	17	20

lowest : 0 1 2 3 4 , highest: 30 31 33 34 46

co: Comborbidities

n	missing	distinct	Info	Mean	Gmd
26082	85	3	0.816	1.348	0.8082

 Value          0     1     2
 Frequency   5035  6931 14116
 Proportion 0.193 0.266 0.541

dcens: Day of censoring for TTR/death days

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
26167	0	29	0.68	27.65	1.254	25	26	28	28	28	29	29

lowest : 1 2 3 4 6 , highest: 29 30 31 32 33

day: Day

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
26167	0	28	0.999	13.88	9.243	2	3	7	14	21	25	27

lowest : 1 2 3 4 5 , highest: 24 25 26 27 28

y: Outcome Status

n	missing	distinct	Info	Mean	Gmd
26167	0	9	0.916	3.773	2.545

lowest : 1.0 1.5 2.0 3.0 4.0 , highest: 4.0 5.0 6.0 7.0 8.0

 Value        1.0   1.5   2.0   3.0   4.0   5.0   6.0   7.0   8.0
 Frequency    930 10765   536   348  2107  3785  1813  5750   133
 Proportion 0.036 0.411 0.020 0.013 0.081 0.145 0.069 0.220 0.005

yprev: Previous Outcome Status

n	missing	distinct	Info	Mean	Gmd
26167	0	8	0.92	3.875	2.529

lowest : 1.0 1.5 2.0 3.0 4.0 , highest: 3.0 4.0 5.0 6.0 7.0

 Value        1.0   1.5   2.0   3.0   4.0   5.0   6.0   7.0
 Frequency    664 10439   415   336  2207  4170  1987  5949
 Proportion 0.025 0.399 0.016 0.013 0.084 0.159 0.076 0.227

gap: Days Since Previous Record

n	missing	distinct	Info	Mean	Gmd
26167	0	8	0.002	1.002	0.003591

lowest : 1 2 3 4 5 , highest: 4 5 6 7 8

 Value          1     2     3     4     5     6     7     8
 Frequency  26153     3     2     2     4     1     1     1
 Proportion 0.999 0.000 0.000 0.000 0.000 0.000 0.000 0.000

seq: Sequential Record Number

n	missing	distinct	Info	Mean	Gmd	.05	.10	.25	.50	.75	.90	.95
26167	0	28	0.999	13.86	9.231	2	3	7	14	21	25	27

lowest : 1 2 3 4 5 , highest: 24 25 26 27 28

Tabulate frequencies of highest outcome category a patient occupies.

w <- d[, .(maxy=max(y, na.rm=TRUE)), by=id]
w[, table(maxy)]

maxy
  4   5   6   7   8 
 89 339 164 326 133

Data Collection Intensity While in Hospital

For patients not dying, the frequency distribution of hospital length of stay is shown below. The 99 is for patients still in hospital at the last follow-up day.

u[is.na(ddeath), table(dlasthospf)]

dlasthospf
  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19 
  5  15  40  51  69  68  45  54  38  38  33  26  22  19  32  25  18  15  12  16 
 20  21  22  23  24  25  26  27  28  29  30  31  32 
 11  28  14   8   9  14  18  13 107  24   2  24   1

Frequencies of occurrences of measurements by day and by hospital length of stay are show in the histograms below. Patients who died are excluded.

w <- d[day <= dlasthospf & is.na(ddeath), ]
ggplot(w, aes(x=day)) + geom_histogram() + facet_wrap(~ dlasthosp)

Summary of Successive State Transitions

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

Category Proportions Over Time

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

dcf$y <- factor(dcf$y, c(1, 1.5, 2 : 8), ordlabels)
gg <- propsPO(y ~ day + tx, nrow=1, data=dcf) +
        theme(legend.position='bottom',
                      axis.text.x=element_text(angle=90, hjust=1))
ggp(gg)

Repeat, dropping any records occurring after the patient arrives at home the first time.

gg <- propsPO(y ~ day + tx, nrow=1, data=dcf[day < pmin(28, dlasthospf),]) +
        theme(legend.position='bottom',
                      axis.text.x=element_text(angle=90, hjust=1))
ggp(gg)

Repeat, dropping any records occurring after first record at home or after time to recovery.

gg <- propsPO(y ~ day + tx, nrow=1,
        data=dcf[day < pmin(28, dlasthospf + 1, ttrecov0 + 1, na.rm=TRUE),]) +
        theme(legend.position='bottom',
                      axis.text.x=element_text(angle=90, hjust=1))
ggp(gg)

Single-Day Analysis

To provide a result similar to the ACTT-1’s pre-specified secondary analysis, we use a proportional odds model for day 15 status in isolation.

w <- d[seq == 1, ]
dd <- datadist(w); options(datadist='dd')

w$y0 <- as.factor(w$y0)
f <- lrm(y15 ~ tx + y0, data=w)
f

Logistic Regression Model

 lrm(formula = y15 ~ tx + y0, data = w)

Frequencies of Responses

   1   2   3   4   5   6   7   8 
 272 219  22  71 118  52 205  92

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 1051	LR χ² 301.96	R² 0.256	C 0.700
max \|∂log L/∂β\| 5×10^-7	d.f. 4	g 1.144	D_xy 0.399
	Pr(>χ²) <0.0001	g_r 3.139	γ 0.462
		g_p 0.252	τ_a 0.329
		Brier 0.179

	β	S.E.	Wald Z	Pr(>\|Z\|)
y≥2	0.2046	0.1684	1.22	0.2243
y≥3	-0.8997	0.1705	-5.28	<0.0001
y≥4	-1.0091	0.1712	-5.89	<0.0001
y≥5	-1.3750	0.1743	-7.89	<0.0001
y≥6	-2.0082	0.1805	-11.12	<0.0001
y≥7	-2.3077	0.1836	-12.57	<0.0001
y≥8	-3.9200	0.2076	-18.88	<0.0001
tx=Remdesivir	-0.3625	0.1107	-3.28	0.0011
y0=5	0.5920	0.1776	3.33	0.0009
y0=6	1.8108	0.2085	8.69	<0.0001
y0=7	2.6637	0.1990	13.39	<0.0001

summary(f)

	Low	High	Effect	S.E.	Lower 0.95	Upper 0.95
Effects Response: `y15`
tx — Placebo:Remdesivir	2	1	0.3625	0.1107	0.1456	0.5794
Odds Ratio	2	1	1.4370		1.1570	1.7850
y0 — 4:5	2	1	-0.5920	0.1776	-0.9401	-0.2439
Odds Ratio	2	1	0.5532		0.3906	0.7836
y0 — 6:5	2	3	1.2190	0.1616	0.9020	1.5360
Odds Ratio	2	3	3.3830		2.4650	4.6440
y0 — 7:5	2	4	2.0720	0.1477	1.7820	2.3610
Odds Ratio	2	4	7.9390		5.9430	10.6000

Time to Recovery Analysis with Cox Model

S <- with(w, Surv(ifelse(! is.na(ttrecov0), ttrecov0, ttrecov1), 
                  ifelse(! is.na(ttrecov0), 1, 0)))
f <- cph(S ~ tx, data=w, x=TRUE, y=TRUE)
f

Cox Proportional Hazards Model

 cph(formula = S ~ tx, data = w, x = TRUE, y = TRUE)

	Model Tests	Discrimination Indexes
Obs 1051	LR χ² 13.73	R² 0.013
Events 751	d.f. 1	D_xy 0.084
Center 0.1376	Pr(>χ²) 0.0002	g 0.136
	Score χ² 13.78	g_r 1.145
	Pr(>χ²) 0.0002

	β	S.E.	Wald Z	Pr(>\|Z\|)
tx=Remdesivir	0.2713	0.0733	3.70	0.0002

summary(f, tx='Placebo')

	Low	High	Effect	S.E.	Lower 0.95	Upper 0.95
Effects Response: `S`
tx — Remdesivir:Placebo	1	2	0.2713	0.0733	0.1276	0.4149
Hazard Ratio	1	2	1.3120		1.1360	1.5140

z <- cox.zph(f, transform='identity')
z

   chisq df     p

tx 5.55 1 0.018 GLOBAL 5.55 1 0.018

plot(z, var='tx')

There is evidence that the treatment effect, when not adjusted for covariates, is not in proportional hazards. The treatment effect on the hazard ratio scale seems to reverse by 20 days as judged by the smoothed scaled Schoenfeld residual plot above. The hazard ratio and 0.95 CL is close to the published value.

Now repeat the model fit adjusting for age and initial state.

f <- cph(S ~ tx + rcs(age, 4) + y0, data=w, x=TRUE, y=TRUE)
f

Cox Proportional Hazards Model

 cph(formula = S ~ tx + rcs(age, 4) + y0, data = w, x = TRUE, 
     y = TRUE)

	Model Tests	Discrimination Indexes
Obs 1051	LR χ² 380.09	R² 0.304
Events 751	d.f. 7	D_xy 0.467
Center -0.9651	Pr(>χ²) 0.0000	g 0.896
	Score χ² 402.58	g_r 2.449
	Pr(>χ²) 0.0000

	β	S.E.	Wald Z	Pr(>\|Z\|)
tx=Remdesivir	0.2641	0.0737	3.58	0.0003
age	0.0039	0.0094	0.42	0.6742
age’	-0.0302	0.0225	-1.34	0.1788
age’’	0.0271	0.1174	0.23	0.8173
y0=5	-0.4829	0.1027	-4.70	<0.0001
y0=6	-1.3062	0.1293	-10.10	<0.0001
y0=7	-1.8941	0.1247	-15.19	<0.0001

anova(f)

	χ²	d.f.	P
Wald Statistics for `S`
tx	12.83	1	0.0003
age	71.11	3	<0.0001
Nonlinear	17.42	2	0.0002
y0	297.28	3	<0.0001
TOTAL	360.17	7	<0.0001

plot(Predict(f, age))
summary(f, tx='Placebo')

	Low	High	Δ	Effect	S.E.	Lower 0.95	Upper 0.95
Effects Response: `S`
age	49	70	21	-0.4637	0.10220	-0.6640	-0.2633
Hazard Ratio	49	70	21	0.6290		0.5148	0.7685
tx — Remdesivir:Placebo	1	2		0.2641	0.07374	0.1196	0.4086
Hazard Ratio	1	2		1.3020		1.1270	1.5050
y0 — 4:5	2	1		0.4829	0.10270	0.2816	0.6842
Hazard Ratio	2	1		1.6210		1.3250	1.9820
y0 — 6:5	2	3		-0.8233	0.10720	-1.0330	-0.6133
Hazard Ratio	2	3		0.4390		0.3558	0.5416
y0 — 7:5	2	4		-1.4110	0.10210	-1.6110	-1.2110
Hazard Ratio	2	4		0.2438		0.1996	0.2979

z <- cox.zph(f, transform='identity')
z

         chisq df      p

tx 2.41 1 0.120 rcs(age, 4) 8.50 3 0.037 y0 94.81 3 <2e-16 GLOBAL 106.72 7 <2e-16

plot(z)

The age effect is extremely nonlinear, with a strong acceleration of effect at age 60. There is less evidence against proportional hazards for treatment, slight non-PH for age, and exceptionally strong non-PH with regard to initial state. Let’s repeat the analysis stratifying on initial state so as to not assume PH for it. We also check for interaction of treatment and initial state.

f <- cph(S ~ rcs(age, 4) + strat(y0) * tx, data=w, x=TRUE, y=TRUE)
anova(f)

	χ²	d.f.	P
Wald Statistics for `S`
age	65.05	3	<0.0001
Nonlinear	15.31	2	0.0005
tx (Factor+Higher Order Factors)	14.55	4	0.0057
All Interactions	4.15	3	0.2453
y0 × tx (Factor+Higher Order Factors)	4.15	3	0.2453
TOTAL NONLINEAR + INTERACTION	19.40	5	0.0016
TOTAL	80.92	7	<0.0001

The treatment interaction is not compelling so let’s drop it.

f <- cph(S ~ rcs(age, 4) + strat(y0) + tx, data=w, x=TRUE, y=TRUE)
summary(f, tx='Placebo')

	Low	High	Δ	Effect	S.E.	Lower 0.95	Upper 0.95
Effects Response: `S`
age	49	70	21	-0.4468	0.10250	-0.64760	-0.2460
Hazard Ratio	49	70	21	0.6397		0.52330	0.7819
tx — Remdesivir:Placebo	1	2		0.2390	0.07379	0.09436	0.3836
Hazard Ratio	1	2		1.2700		1.09900	1.4680

anova(f)

	χ²	d.f.	P
Wald Statistics for `S`
age	65.81	3	<0.0001
Nonlinear	15.27	2	0.0005
tx	10.49	1	0.0012
TOTAL	76.90	4	<0.0001

z <- cox.zph(f, transform='identity')
z

        chisq df     p

rcs(age, 4) 10.06 3 0.018 tx 3.41 1 0.065 GLOBAL 12.96 4 0.011

plot(z, var='tx')

The treatment effect remains in non-proportional hazard, and there is also non-PH for age. There is an apparent reversal of treatment effect at 15d.

Take a look at the time-varying effect of initial state to better understand its non-PH. First get transformed stratified K-M curves ignoring all other variables, then fit initial state in a Cox model but estimate the time-varying log hazard ratio.

f <- npsurv(S ~ y0, data=w)
survplot(f, fun=function(y) log(-log(y)))
f <- cph(S ~ y0, data=w, x=TRUE, y=TRUE)
z <- cox.zph(f, transform='identity', terms=FALSE)
z

         chisq df       p
y0=5    33.397  1 7.5e-09
y0=6     0.319  1    0.57
y0=7    97.703  1 < 2e-16
GLOBAL 101.810  3 < 2e-16

par(mfrow=c(2,2))
plot(z)

The most non-parallelism comes from initial state 7 (invasive ventilation), which wears off over time.

Serial Conditional Discrete Time Markov Ordinal Model

Here we use a first order discrete time 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. The 1-day lagged outcome status is variable yprev. Note that the gaps between measurement occurs largely once the patient is home. Follow-up is curtailed at 30d.

Start by examining frequencies of covariate combinations.

a <- d
a[, table(yprev, y)]

     y
yprev    1  1.5    2    3    4    5    6    7    8
  1      3  658    1    0    2    0    0    0    0
  1.5  912 8990  524    0    2    3    0    0    8
  2      1  395    0    0   11    6    2    0    0
  3      2   46    1  260   18    8    1    0    0
  4      6  375    7   62 1601  137   15    1    3
  5      5  262    3   23  415 3215  170   57   20
  6      1   30    0    3   52  294 1439  145   23
  7      0    9    0    0    6  122  186 5547   79

a[, ggfreqScatter(yprev, y)]
a[, ggfreqScatter(gap, yprev)]

Fit a fairly saturated ordinary proportional odds 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
covariate adjustment for a flexible function of age (in addition to baseline state)

This model does not use interval censoring for y=1.5 but treats this category as one between y=1 and y=1.

a$yprev <- factor(a$yprev)
dd <- datadist(a); options(datadist='dd')
f <- lrm(y ~ yprev * tx + rcs(day, 4) * (tx + yprev) + rcs(age, 5),
          data=a)
anova(f)

	χ²	d.f.	P
Wald Statistics for `y`
yprev (Factor+Higher Order Factors)	15637.37	35	<0.0001
All Interactions	136.23	28	<0.0001
tx (Factor+Higher Order Factors)	13.70	11	0.2498
All Interactions	7.31	10	0.6958
day (Factor+Higher Order Factors)	340.41	27	<0.0001
All Interactions	139.91	24	<0.0001
Nonlinear (Factor+Higher Order Factors)	200.45	18	<0.0001
age	103.33	4	<0.0001
Nonlinear	7.35	3	0.0617
yprev × tx (Factor+Higher Order Factors)	1.49	7	0.9827
tx × day (Factor+Higher Order Factors)	6.19	3	0.1025
Nonlinear	0.90	2	0.6369
Nonlinear Interaction : f(A,B) vs. AB	0.90	2	0.6369
yprev × day (Factor+Higher Order Factors)	135.49	21	<0.0001
Nonlinear	80.18	14	<0.0001
Nonlinear Interaction : f(A,B) vs. AB	80.18	14	<0.0001
TOTAL NONLINEAR	206.77	21	<0.0001
TOTAL INTERACTION	140.93	31	<0.0001
TOTAL NONLINEAR + INTERACTION	233.59	36	<0.0001
TOTAL	16074.66	46	<0.0001

There is little evidence for the treatment effect depending on the previous state. Remove that interaction and re-fit. Plot how the effect of previous state and treatment depends on time.

f <- lrm(y ~ rcs(day, 4) * (tx + yprev) + rcs(age, 5), data=a)
f

Logistic Regression Model

 lrm(formula = y ~ rcs(day, 4) * (tx + yprev) + rcs(age, 5), data = a)

Frequencies of Responses

     1   1.5     2     3     4     5     6     7     8 
   930 10765   536   348  2107  3785  1813  5750   133

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 26167	LR χ² 51497.39	R² 0.893	C 0.952
max \|∂log L/∂β\| 7×10^-11	d.f. 39	g 7.042	D_xy 0.904
	Pr(>χ²) <0.0001	g_r 1143.542	γ 0.914
		g_p 0.487	τ_a 0.677
		Brier 0.030

	β	S.E.	Wald Z	Pr(>\|Z\|)
y≥1.5	2.8969	2.6211	1.11	0.2691
y≥2	-2.6051	2.6210	-0.99	0.3203
y≥3	-3.1880	2.6211	-1.22	0.2239
y≥4	-3.6700	2.6212	-1.40	0.1615
y≥5	-6.0997	2.6217	-2.33	0.0200
y≥6	-9.9883	2.6224	-3.81	0.0001
y≥7	-12.6845	2.6231	-4.84	<0.0001
y≥8	-19.2076	2.6251	-7.32	<0.0001
day	0.0052	0.3084	0.02	0.9866
day’	-0.1078	0.7099	-0.15	0.8793
day’’	0.3382	1.7969	0.19	0.8507
tx=Remdesivir	-0.2077	0.1033	-2.01	0.0443
yprev=1.5	0.3170	2.6358	0.12	0.9043
yprev=2	3.0293	4.4233	0.68	0.4934
yprev=3	3.2314	2.7107	1.19	0.2332
yprev=4	5.1244	2.6138	1.96	0.0499
yprev=5	8.3390	2.6131	3.19	0.0014
yprev=6	11.2879	2.6146	4.32	<0.0001
yprev=7	15.6574	2.6177	5.98	<0.0001
age	0.0003	0.0055	0.06	0.9543
age’	0.0238	0.0236	1.01	0.3148
age’’	-0.0396	0.1316	-0.30	0.7632
age’’’	-0.0572	0.2494	-0.23	0.8185
day × tx=Remdesivir	0.0082	0.0187	0.44	0.6591
day’ × tx=Remdesivir	0.0218	0.0648	0.34	0.7365
day’’ × tx=Remdesivir	-0.0910	0.1877	-0.48	0.6280
day × yprev=1.5	-0.0793	0.3117	-0.25	0.7993
day’ × yprev=1.5	0.1898	0.7176	0.26	0.7914
day’’ × yprev=1.5	-0.5232	1.8145	-0.29	0.7731
day × yprev=2	0.0984	0.5278	0.19	0.8522
day’ × yprev=2	-1.4550	1.1893	-1.22	0.2212
day’’ × yprev=2	5.0168	2.9453	1.70	0.0885
day × yprev=3	-0.0046	0.3233	-0.01	0.9887
day’ × yprev=3	0.0536	0.7530	0.07	0.9433
day’’ × yprev=3	-0.1502	1.9124	-0.08	0.9374
day × yprev=4	-0.1713	0.3089	-0.55	0.5793
day’ × yprev=4	0.6241	0.7132	0.88	0.3815
day’’ × yprev=4	-1.7133	1.8089	-0.95	0.3436
day × yprev=5	-0.2074	0.3087	-0.67	0.5016
day’ × yprev=5	0.6383	0.7124	0.90	0.3703
day’’ × yprev=5	-1.6423	1.8068	-0.91	0.3634
day × yprev=6	-0.0775	0.3089	-0.25	0.8018
day’ × yprev=6	0.1611	0.7144	0.23	0.8216
day’’ × yprev=6	-0.2824	1.8153	-0.16	0.8764
day × yprev=7	-0.0585	0.3093	-0.19	0.8500
day’ × yprev=7	0.1676	0.7152	0.23	0.8148
day’’ × yprev=7	-0.4150	1.8160	-0.23	0.8192

anova(f)

	χ²	d.f.	P
Wald Statistics for `y`
day (Factor+Higher Order Factors)	339.67	27	<0.0001
All Interactions	139.47	24	<0.0001
Nonlinear (Factor+Higher Order Factors)	200.10	18	<0.0001
tx (Factor+Higher Order Factors)	12.22	4	0.0158
All Interactions	5.83	3	0.1203
yprev (Factor+Higher Order Factors)	15639.11	28	<0.0001
All Interactions	134.77	21	<0.0001
age	103.37	4	<0.0001
Nonlinear	7.22	3	0.0653
day × tx (Factor+Higher Order Factors)	5.83	3	0.1203
Nonlinear	0.88	2	0.6451
Nonlinear Interaction : f(A,B) vs. AB	0.88	2	0.6451
day × yprev (Factor+Higher Order Factors)	134.77	21	<0.0001
Nonlinear	80.09	14	<0.0001
Nonlinear Interaction : f(A,B) vs. AB	80.09	14	<0.0001
TOTAL NONLINEAR	206.28	21	<0.0001
TOTAL INTERACTION	139.47	24	<0.0001
TOTAL NONLINEAR + INTERACTION	232.13	29	<0.0001
TOTAL	16074.50	39	<0.0001

ggplot(Predict(f, day, tx))
ggplot(Predict(f, day, yprev))

There is a strange time-varying effect of being previously in state 2. Let’s examine the frequency distribution by time. See that there’s not many previous states of 2 except for days 14-24.

a[, table(day, yprev)]

    yprev
day    1 1.5   2   3   4   5   6   7
  1    0   0   0   0 138 435 193 285
  2    0   0   0   1 123 405 199 317
  3    2  12   0   3 129 376 172 339
  4    0  51   0   3 127 332 157 352
  5    1  97   0   9 137 287 137 347
  6    0 166   2   7 120 261 116 337
  7    1 233   0  12 113 213 106 312
  8    1 277   1  14 109 178  93 306
  9    1 330   2  16  92 152  89 291
  10   3 362   1  18  87 143  70 283
  11   1 400   5  23  73 127  68 268
  12   3 430   1  20  67 114  65 256
  13  10 437  11  20  69 102  65 238
  14  34 422  22  21  70 101  57 220
  15 115 299  77  22  68 102  47 201
  16  64 409  48  18  67  90  40 187
  17  37 482  26  16  55  90  41 171
  18  11 542  10  14  60  83  36 155
  19   6 562   9  12  54  85  32 147
  20   9 568  12  10  64  73  31 136
  21  39 542  22  11  52  76  31 125
  22 165 375  67  13  59  67  24 124
  23  84 513  40  11  54  52  23 112
  24  40 568  41   9  52  48  20 105
  25  23 614  14   9  52  50  22  94
  26  13 620   2  11  40  45  21  91
  27   1 600   2   7  44  41  15  82
  28   0 528   0   6  32  42  17  68

See if there is any evidence for a time-dependent previous state effect other than for previous state 2.

w <- a[yprev != '2', ]
w[, yprev := yprev[drop=TRUE]]
g <- lrm(y ~ rcs(day, 4) * (tx + yprev) + rcs(age, 5), data=w)

anova(g)

	χ²	d.f.	P
Wald Statistics for `y`
day (Factor+Higher Order Factors)	306.53	24	<0.0001
All Interactions	100.20	21	<0.0001
Nonlinear (Factor+Higher Order Factors)	172.23	16	<0.0001
tx (Factor+Higher Order Factors)	12.28	4	0.0154
All Interactions	5.74	3	0.1250
yprev (Factor+Higher Order Factors)	15520.00	24	<0.0001
All Interactions	95.48	18	<0.0001
age	103.06	4	<0.0001
Nonlinear	7.17	3	0.0667
day × tx (Factor+Higher Order Factors)	5.74	3	0.1250
Nonlinear	0.75	2	0.6889
Nonlinear Interaction : f(A,B) vs. AB	0.75	2	0.6889
day × yprev (Factor+Higher Order Factors)	95.48	18	<0.0001
Nonlinear	46.72	12	<0.0001
Nonlinear Interaction : f(A,B) vs. AB	46.72	12	<0.0001
TOTAL NONLINEAR	178.46	19	<0.0001
TOTAL INTERACTION	100.20	21	<0.0001
TOTAL NONLINEAR + INTERACTION	201.85	26	<0.0001
TOTAL	15936.33	35	<0.0001

ggplot(Predict(g, day, yprev))

Though lessened, there is still some evidence for time by previous state interaction. The nonlinear previous state by time interaction causes an instability that will ultimately distort state occupancy probabilities at later times. With the amount of explained outcome variation due to this interaction being 87/15237 of the total effect of previous state, re-fit the model excluding this interaction. Use the new model to get daily treatment contrasts

f <- lrm(y ~ yprev + rcs(day, 4) * tx + rcs(age, 5), data=a)
f

Logistic Regression Model

 lrm(formula = y ~ yprev + rcs(day, 4) * tx + rcs(age, 5), data = a)

Frequencies of Responses

     1   1.5     2     3     4     5     6     7     8 
   930 10765   536   348  2107  3785  1813  5750   133

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 26167	LR χ² 51367.34	R² 0.893	C 0.952
max \|∂log L/∂β\| 7×10^-11	d.f. 18	g 6.991	D_xy 0.903
	Pr(>χ²) <0.0001	g_r 1086.958	γ 0.913
		g_p 0.486	τ_a 0.676
		Brier 0.031

	β	S.E.	Wald Z	Pr(>\|Z\|)
y≥1.5	3.5539	0.2543	13.98	<0.0001
y≥2	-1.9109	0.2543	-7.51	<0.0001
y≥3	-2.4912	0.2555	-9.75	<0.0001
y≥4	-2.9697	0.2566	-11.57	<0.0001
y≥5	-5.3915	0.2607	-20.68	<0.0001
y≥6	-9.2260	0.2659	-34.70	<0.0001
y≥7	-11.9289	0.2718	-43.89	<0.0001
y≥8	-18.4854	0.2888	-64.00	<0.0001
yprev=1.5	-0.3220	0.1189	-2.71	0.0068
yprev=2	0.3187	0.1843	1.73	0.0837
yprev=3	3.3764	0.1473	22.92	<0.0001
yprev=4	4.3627	0.1309	33.34	<0.0001
yprev=5	7.2692	0.1394	52.14	<0.0001
yprev=6	10.7192	0.1544	69.42	<0.0001
yprev=7	15.2709	0.1671	91.38	<0.0001
day	-0.1150	0.0136	-8.44	<0.0001
day’	0.2561	0.0469	5.46	<0.0001
day’’	-0.6344	0.1357	-4.68	<0.0001
tx=Remdesivir	-0.1910	0.1030	-1.86	0.0636
age	-0.0003	0.0055	-0.05	0.9629
age’	0.0276	0.0236	1.17	0.2417
age’’	-0.0691	0.1312	-0.53	0.5983
age’’’	0.0044	0.2487	0.02	0.9858
day × tx=Remdesivir	0.0068	0.0186	0.37	0.7124
day’ × tx=Remdesivir	0.0256	0.0645	0.40	0.6919
day’’ × tx=Remdesivir	-0.1079	0.1869	-0.58	0.5636

anova(f)

	χ²	d.f.	P
Wald Statistics for `y`
yprev	15698.66	7	<0.0001
day (Factor+Higher Order Factors)	207.77	6	<0.0001
All Interactions	4.70	3	0.1955
Nonlinear (Factor+Higher Order Factors)	90.30	4	<0.0001
tx (Factor+Higher Order Factors)	10.97	4	0.0269
All Interactions	4.70	3	0.1955
age	96.57	4	<0.0001
Nonlinear	6.31	3	0.0973
day × tx (Factor+Higher Order Factors)	4.70	3	0.1955
Nonlinear	1.29	2	0.5235
Nonlinear Interaction : f(A,B) vs. AB	1.29	2	0.5235
TOTAL NONLINEAR	96.06	7	<0.0001
TOTAL NONLINEAR + INTERACTION	99.45	8	<0.0001
TOTAL	16136.02	18	<0.0001

ggplot(Predict(f, day, tx))

contrast(f, list(tx='Placebo',    day=1:28), 
            list(tx='Remdesivir', day=1:28))

     yprev day age     Contrast       S.E.        Lower     Upper     Z
 1     1.5   1  59  0.184179256 0.08702031  0.013622591 0.3547359  2.12
 2     1.5   2  59  0.177332650 0.07234119  0.035546529 0.3191188  2.45
 3     1.5   3  59  0.170445151 0.05989887  0.053045514 0.2878448  2.85
 4*    1.5   4  59  0.163312292 0.05128259  0.062800259 0.2638243  3.18
 5*    1.5   5  59  0.155688712 0.04779505  0.062012127 0.2493653  3.26
 6*    1.5   6  59  0.147329052 0.04907197  0.051149765 0.2435083  3.00
 7*    1.5   7  59  0.137987951 0.05299716  0.034115435 0.2418605  2.60
 8*    1.5   8  59  0.127420050 0.05719942  0.015311241 0.2395289  2.23
 9*    1.5   9  59  0.115379987 0.05995859 -0.002136692 0.2328967  1.92
10*    1.5  10  59  0.101622404 0.06025491 -0.016475059 0.2197199  1.69
11     1.5  11  59  0.086074595 0.05794018 -0.027486067 0.1996353  1.49
12*    1.5  12  59  0.069354478 0.05432593 -0.037122394 0.1758314  1.28
13*    1.5  13  59  0.052252628 0.05127860 -0.048251591 0.1527568  1.02
14*    1.5  14  59  0.035559617 0.05049058 -0.063400096 0.1345193  0.70
15*    1.5  15  59  0.020066020 0.05260285 -0.083033678 0.1231657  0.38
16*    1.5  16  59  0.006562410 0.05677327 -0.104711159 0.1178360  0.12
17*    1.5  17  59 -0.004160640 0.06128563 -0.124278267 0.1159570 -0.07
18*    1.5  18  59 -0.011312555 0.06436757 -0.137470674 0.1148456 -0.18
19*    1.5  19  59 -0.014315298 0.06481113 -0.141342778 0.1127122 -0.22
20*    1.5  20  59 -0.013440966 0.06280705 -0.136540531 0.1096586 -0.21
21*    1.5  21  59 -0.009174195 0.05935644 -0.125510678 0.1071623 -0.15
22*    1.5  22  59 -0.001999619 0.05608535 -0.111924884 0.1079256 -0.04
23*    1.5  23  59  0.007598128 0.05518332 -0.100559190 0.1157554  0.14
24*    1.5  24  59  0.019134412 0.05870453 -0.095924346 0.1341932  0.33
25*    1.5  25  59  0.032124598 0.06735515 -0.099889078 0.1641383  0.48
26*    1.5  26  59  0.046084054 0.08029832 -0.111297764 0.2034659  0.57
27*    1.5  27  59  0.060528143 0.09611867 -0.127860998 0.2489173  0.63
28*    1.5  28  59  0.075053004 0.11358786 -0.147575108 0.2976811  0.66
     Pr(>|z|)
 1     0.0343
 2     0.0142
 3     0.0044
 4*    0.0014
 5*    0.0011
 6*    0.0027
 7*    0.0092
 8*    0.0259
 9*    0.0543
10*    0.0917
11     0.1374
12*    0.2017
13*    0.3082
14*    0.4813
15*    0.7029
16*    0.9080
17*    0.9459
18*    0.8605
19*    0.8252
20*    0.8305
21*    0.8772
22*    0.9716
23*    0.8905
24*    0.7445
25*    0.6334
26*    0.5660
27*    0.5289
28*    0.5088

Redundant contrasts are denoted by *

Confidence intervals are 0.95 individual intervals

The treatment effect, in terms of state transition probabilities, is concentrated in the first 11 days. This is in line with what was seen in the earlier time-dependent treatment log hazard ratio plot from the time to recovery analysis.

Proportional Odds Assumption

Let’s check the PO assumption with respect to treatment, time, and yprev. We do this by fiting a binary logistic model for all Y-cutoffs and tracking regression coefficients against cutoff. For simplicity we force linearity remove interactions.

n <- c('yprev=1.5', 'yprev=2', 'yprev=3', 'yprev=4', 'yprev=5', 'yprev=6', 'yprev=7', 'day', 'tx=Remdesivir')
z <- se <- matrix(NA, nrow=8, ncol=length(n))
cuts <- c(1.5, 2:8)
dimnames(z) <- dimnames(se) <- list(paste0('y>=', cuts), n)
i <- 0
for(ycut in cuts) {
  i <- i + 1
    f <- lrm(y >= ycut ~ yprev + day + tx, data=a, tol=1e-14, maxit=50)
    z [i, n] <- coef(f)[n]
    se[i, n] <- sqrt(diag(vcov(f)[n, n])) 
}
round(z, 2)

       yprev=1.5 yprev=2 yprev=3 yprev=4 yprev=5 yprev=6 yprev=7   day
y>=1.5     -3.08    0.59   -0.53    0.05    0.66    1.50    9.35 -0.08
y>=2        2.49    2.36    7.22    7.02    8.15    9.61   11.94  0.01
y>=3       -0.90    2.76    7.51    7.24    8.33    9.80   12.18 -0.01
y>=4       -0.91    2.75    3.29    7.03    8.20    9.67   12.15 -0.02
y>=5        6.32    9.27    9.46   10.34   14.46   15.97   18.93 -0.04
y>=6        5.20    7.08    6.43    7.32    9.20   13.46   15.88 -0.04
y>=7        6.43   -0.02   -0.15    7.05    9.30   10.89   16.25 -0.03
y>=8        6.12    0.01    0.08    6.82    8.13    9.02    9.11  0.02
       tx=Remdesivir
y>=1.5         -0.05
y>=2           -0.08
y>=3           -0.06
y>=4           -0.07
y>=5           -0.15
y>=6           -0.18
y>=7           -0.18
y>=8           -0.10

round(se, 2)

       yprev=1.5 yprev=2 yprev=3 yprev=4 yprev=5 yprev=6 yprev=7  day
y>=1.5      0.58    1.16    0.92    0.71    0.73    1.16   26.78 0.01
y>=2        0.58    0.62    0.60    0.58    0.58    0.61    0.67 0.00
y>=3        0.76    0.75    0.73    0.71    0.71    0.73    0.78 0.01
y>=4        0.76    0.75    0.74    0.71    0.71    0.73    0.78 0.00
y>=5       28.85   28.85   28.85   28.85   28.85   28.85   28.85 0.00
y>=6       19.52   19.53   19.54   19.52   19.52   19.52   19.52 0.01
y>=7       35.67   57.51   61.48   35.67   35.67   35.67   35.67 0.01
y>=8       29.70   47.88   51.15   29.70   29.70   29.70   29.70 0.01
       tx=Remdesivir
y>=1.5          0.07
y>=2            0.06
y>=3            0.08
y>=4            0.07
y>=5            0.07
y>=6            0.08
y>=7            0.09
y>=8            0.18

Keeping in mind that many of the other entries in the first table presented above have high standard errors (shown in the second table), the variable that is most systematically violating PO is time, just as with VIOLET and ORCHID. There is some indication that treatment has a larger effect on outcome levels 5-7.

Let’s fit a partial proportional odds model to handle time and also to assess PO for treatment. Drop the yprev by day interaction and nonlinear terms for now.

require(VGAM)

f <- vgam(ordered(y) ~ yprev + day + tx + age,
          cumulative(parallel=TRUE, reverse=TRUE),
          data=a)
h <- vgam(ordered(y) ~ yprev + day + tx + age,
          cumulative(parallel=FALSE ~ tx + day, reverse=TRUE),
          data=a)

setdiff(names(coef(h)), names(coef(f)))

 [1] "day:1"          "day:2"          "day:3"          "day:4"         
 [5] "day:5"          "day:6"          "day:7"          "day:8"         
 [9] "txRemdesivir:1" "txRemdesivir:2" "txRemdesivir:3" "txRemdesivir:4"
[13] "txRemdesivir:5" "txRemdesivir:6" "txRemdesivir:7" "txRemdesivir:8"

setdiff(names(coef(f)), names(coef(h)))

[1] "day"          "txRemdesivir"

cbind(beta=round(coef(h),3), SE=round(sqrt(diag(vcov(h))), 3))

                 beta    SE
(Intercept):1   4.155 0.176
(Intercept):2  -2.500 0.121
(Intercept):3  -2.484 0.121
(Intercept):4  -2.491 0.121
(Intercept):5  -2.990 0.124
(Intercept):6  -4.326 0.125
(Intercept):7  -5.142 0.126
(Intercept):8  -9.941 0.242
yprev1.5       -0.357 0.122
yprev2          0.325 0.188
yprev3          2.481 0.121
yprev4          2.678 0.121
yprev5          3.482 0.121
yprev6          4.256 0.121
yprev7          5.273 0.121
day:1          -0.078 0.006
day:2          -0.019 0.000
day:3          -0.027 0.000
day:4          -0.034 0.001
day:5          -0.038 0.002
day:6          -0.014 0.002
day:7           0.004 0.002
day:8           0.018 0.013
txRemdesivir:1 -0.036 0.069
txRemdesivir:2  0.018 0.000
txRemdesivir:3  0.015 0.000
txRemdesivir:4  0.004 0.010
txRemdesivir:5 -0.106 0.026
txRemdesivir:6 -0.177 0.033
txRemdesivir:7 -0.188 0.036
txRemdesivir:8 -0.063 0.205
age             0.003 0.000

Fit a Bayesian constrained partial proportional odds model with time and age as linear effects and allowing non-PO for time and treatment. The constraint is that the departure from PO is linear in y. Use proper interval censoring for patients known to be at home for which it is not known whether y=1 vs. y=2. For now force the time effect to be linear.

stanSet()

a[, yl := ifelse(y == 1.5, 1, y)]
a[, yu := ifelse(y == 1.5, 2, y)]
bppo <- blrm(Ocens(yl, yu) ~ yprev + day * tx + age,
             ~ tx + day, cppo = function(y) y,
          data=a, file='bppo.rds', refresh=100)
bppo

Bayesian Constrained Partial Proportional Odds Ordinal Logistic Model

Dirichlet Priors With Concentration Parameter 0.278 for Intercepts

 blrm(formula = Ocens(yl, yu) ~ yprev + day * tx + age, ppo = ~tx + 
     day, cppo = function(y) y, data = a, refresh = 100, file = "bppo.rds")

Frequencies of Responses

    1    2    3    4    5    6    7    8 
  930  536  348 2107 3785 1813 5750  133

	Mixed Calibration/ Discrimination Indexes	Discrimination Indexes	Rank Discrim. Indexes
Obs 26167	B 0.034 [0.033, 0.034]	g 8.036 [7.834, 8.234]	C 0.964 [0.963, 0.964]
Censored 10765		g_p 0.482 [0.481, 0.483]	D_xy 0.927 [0.925, 0.929]
L=10765		EV 0.841 [0.835, 0.847]
Draws 4000		v 52.772 [49.969, 55.067]
Chains 4		vp 0.208 [0.206, 0.209]
p 11

	Mode β	Mean β	Median β	S.E.	Lower	Upper	Pr(β>0)	Symmetry
y≥2	-1.1737	-1.3711	-1.3413	0.6313	-2.6245	-0.2097	0.0075	0.84
y≥3	-6.0610	-6.3352	-6.2848	0.6164	-7.6110	-5.2197	0.0000	0.84
y≥4	-6.4389	-6.7135	-6.6640	0.6144	-7.9376	-5.5441	0.0000	0.85
y≥5	-8.7561	-9.0362	-8.9968	0.6137	-10.2338	-7.8393	0.0000	0.87
y≥6	-12.4526	-12.7368	-12.6909	0.6145	-13.9607	-11.5730	0.0000	0.86
y≥7	-15.0427	-15.3312	-15.2921	0.6186	-16.4934	-14.0873	0.0000	0.87
y≥8	-21.5045	-21.8020	-21.7651	0.6290	-23.0072	-20.5591	0.0000	0.87
yprev=1.5	0.0395	0.2205	0.1914	0.6129	-0.9804	1.3582	0.6228	1.20
yprev=2	2.5380	2.7911	2.7599	0.6493	1.5750	4.1265	1.0000	1.18
yprev=3	5.8124	6.0762	6.0326	0.6098	4.9820	7.3274	1.0000	1.20
yprev=4	6.9164	7.1920	7.1469	0.6057	6.0235	8.3610	1.0000	1.22
yprev=5	9.7847	10.0655	10.0173	0.6070	9.0012	11.3501	1.0000	1.20
yprev=6	13.1623	13.4470	13.4059	0.6100	12.3055	14.6472	1.0000	1.17
yprev=7	17.6585	17.9543	17.9159	0.6164	16.7718	19.1359	1.0000	1.17
day	-0.0134	-0.0130	-0.0130	0.0043	-0.0211	-0.0045	0.0013	1.02
tx=Remdesivir	-0.1497	-0.1549	-0.1543	0.0782	-0.3126	-0.0042	0.0260	1.00
age	0.0101	0.0101	0.0101	0.0012	0.0077	0.0125	1.0000	1.00
day × tx=Remdesivir	0.0066	0.0068	0.0067	0.0047	-0.0029	0.0162	0.9255	1.01
tx=Remdesivir x f(y)	-0.0361	-0.0314	-0.0326	0.0516	-0.1327	0.0689	0.2702	1.02
day x f(y)	-0.0189	-0.0195	-0.0195	0.0038	-0.0269	-0.0121	0.0000	1.01

As with VIOLET and ORCHID, the major non-PO effect is with respect to time, with must less evidence for non-PO for treatment. Refit the model allowing only non-PO for (linear) time even though the main time effect is nonlinear. Also adjust for (nonlinear) age.

bppodc <- blrm(Ocens(yl, yu) ~ yprev + rcs(day, 4) * tx + rcs(age, 5),
             ~ day, cppo = function(y) y,
          data=a, file='bppodc.rds', refresh=100)
stanDx(bppodc)

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                   2417 1.000
y>=3                   2191 1.000
y>=4                   2030 1.001
y>=5                   2050 1.000
y>=6                   2435 1.000
y>=7                   2528 1.001
y>=8                   2861 1.001
yprev=1.5              2198 1.000
yprev=2                4745 1.000
yprev=3                3177 1.001
yprev=4                2212 1.002
yprev=5                2135 1.001
yprev=6                2858 1.001
yprev=7                3608 0.999
day                    3672 1.000
day'                   4931 1.000
day''                  5182 0.999
tx=Remdesivir          4468 1.000
age                    5510 1.001
age'                   5580 1.000
age''                  4865 1.000
age'''                 5209 1.000
day * tx=Remdesivir    4742 1.000
day' * tx=Remdesivir   5725 0.999
day'' * tx=Remdesivir  5599 1.000
day x f(y)             3493 1.001

stanDxplot(bppodc)

bppodc

Bayesian Constrained Partial Proportional Odds Ordinal Logistic Model

Dirichlet Priors With Concentration Parameter 0.278 for Intercepts

 blrm(formula = Ocens(yl, yu) ~ yprev + rcs(day, 4) * tx + rcs(age, 
     5), ppo = ~day, cppo = function(y) y, data = a, refresh = 100, 
     file = "bppodc.rds")

Frequencies of Responses

    1    2    3    4    5    6    7    8 
  930  536  348 2107 3785 1813 5750  133

	Mixed Calibration/ Discrimination Indexes	Discrimination Indexes	Rank Discrim. Indexes
Obs 26167	B 0.033 [0.033, 0.033]	g 8.058 [7.887, 8.247]	C 0.966 [0.965, 0.967]
Censored 10765		g_p 0.483 [0.482, 0.483]	D_xy 0.932 [0.929, 0.934]
L=10765		EV 0.842 [0.836, 0.847]
Draws 4000		v 52.992 [50.939, 55.589]
Chains 4		vp 0.208 [0.206, 0.209]
p 18

	Mode β	Mean β	Median β	S.E.	Lower	Upper	Pr(β>0)	Symmetry
y≥2	-0.3256	-0.4592	-0.4300	0.6894	-1.7670	0.9126	0.2602	0.84
y≥3	-5.3321	-5.4592	-5.4221	0.6738	-6.7534	-4.0999	0.0000	0.85
y≥4	-5.7558	-5.8862	-5.8533	0.6709	-7.1595	-4.5295	0.0000	0.86
y≥5	-8.1318	-8.2625	-8.2243	0.6705	-9.5575	-6.9351	0.0000	0.86
y≥6	-11.8988	-12.0268	-11.9874	0.6713	-13.3017	-10.6628	0.0000	0.86
y≥7	-14.5419	-14.6690	-14.6313	0.6730	-16.0453	-13.4099	0.0000	0.88
y≥8	-21.0513	-21.1842	-21.1477	0.6816	-22.5979	-19.9375	0.0000	0.88
yprev=1.5	0.1296	0.2718	0.2366	0.6232	-0.8464	1.5745	0.6577	1.21
yprev=2	2.6786	2.8023	2.7695	0.6620	1.6219	4.1813	1.0000	1.17
yprev=3	6.0374	6.1627	6.1235	0.6210	5.0376	7.4118	1.0000	1.21
yprev=4	7.0791	7.2088	7.1685	0.6179	6.0564	8.4328	1.0000	1.21
yprev=5	9.9557	10.0827	10.0423	0.6186	8.9728	11.3574	1.0000	1.22
yprev=6	13.3842	13.5051	13.4654	0.6221	12.3537	14.7492	1.0000	1.21
yprev=7	17.9329	18.0530	18.0126	0.6279	16.8958	19.3063	1.0000	1.21
day	-0.1089	-0.1088	-0.1091	0.0150	-0.1388	-0.0793	0.0000	0.98
day’	0.2821	0.2797	0.2796	0.0516	0.1753	0.3799	1.0000	1.02
day’’	-0.7027	-0.6947	-0.6940	0.1512	-0.9854	-0.3904	0.0000	0.99
tx=Remdesivir	-0.1642	-0.1625	-0.1631	0.1028	-0.3643	0.0342	0.0565	0.97
age	0.0004	0.0004	0.0004	0.0064	-0.0120	0.0132	0.5308	0.99
age’	0.0232	0.0230	0.0229	0.0272	-0.0313	0.0756	0.8050	0.95
age’’	-0.0478	-0.0475	-0.0481	0.1505	-0.3575	0.2304	0.3820	1.00
age’’’	-0.0158	-0.0150	-0.0114	0.2828	-0.5428	0.5516	0.4822	0.99
day × tx=Remdesivir	-0.0019	-0.0022	-0.0023	0.0196	-0.0415	0.0351	0.4542	1.01
day’ × tx=Remdesivir	0.0614	0.0628	0.0622	0.0707	-0.0752	0.1996	0.8168	1.01
day’’ × tx=Remdesivir	-0.2050	-0.2089	-0.2067	0.2076	-0.6254	0.1798	0.1582	0.97
day x f(y)	-0.0128	-0.0123	-0.0123	0.0041	-0.0204	-0.0045	0.0018	0.98

anova(bppodc)

	REV	Lower	Upper	d.f.
Relative Explained Variation for `Ocens(yl, yu)`. Approximate total model Wald χ² used in denominators of REV:12558.6 [12144.5, 13062.7].
yprev	0.941	0.933	0.948	7
day (Factor+Higher Order Factors)	0.008	0.006	0.012	6
All Interactions	0.000	0.000	0.001	3
Nonlinear (Factor+Higher Order Factors)	0.008	0.005	0.012	4
tx (Factor+Higher Order Factors)	0.001	0.000	0.002	4
All Interactions	0.000	0.000	0.001	3
age	0.006	0.004	0.010	4
Nonlinear	0.000	0.000	0.001	3
day × tx (Factor+Higher Order Factors)	0.000	0.000	0.001	3
Nonlinear	0.000	0.000	0.001	2
Nonlinear Interaction : f(A,B) vs. AB	0.000	0.000	0.001	2
TOTAL NONLINEAR	0.008	0.006	0.012	7
TOTAL NONLINEAR + INTERACTION	0.008	0.006	0.012	8
TOTAL	0.990	0.987	0.994	18

ggplot(Predict(bppodc, day, tx))

# Choice of ycut won't matter for transition ORs for treatment because tx is in PO
k <- contrast(bppodc, ycut=4,
              list(tx='Remdesivir', day=1:28), list(tx='Placebo', day=1:28))
k$cnames <- 'y >= 4'
k

           yprev day age     Contrast       S.E.       Lower        Upper
 1  y >= 4   1.5   1  59 -0.164718262 0.08650997 -0.33782070 -0.002784024
 2  y >= 4   1.5   2  59 -0.166950617 0.07182736 -0.31795872 -0.036617261
 3  y >= 4   1.5   3  59 -0.169082528 0.06003179 -0.28924617 -0.053996492
 4* y >= 4   1.5   4  59 -0.170611776 0.05293808 -0.27719054 -0.069094852
 5* y >= 4   1.5   5  59 -0.170935696 0.05157570 -0.27707360 -0.070057793
 6* y >= 4   1.5   6  59 -0.169451625 0.05481486 -0.28001054 -0.063008998
 7* y >= 4   1.5   7  59 -0.165556900 0.06008602 -0.28416971 -0.050833589
 8* y >= 4   1.5   8  59 -0.158648855 0.06503332 -0.28793354 -0.033674372
 9* y >= 4   1.5   9  59 -0.148124828 0.06807096 -0.28678030 -0.021918483
10* y >= 4   1.5  10  59 -0.133382155 0.06830520 -0.27071702 -0.005791486
11  y >= 4   1.5  11  59 -0.114152438 0.06576022 -0.24579458  0.009393732
12* y >= 4   1.5  12  59 -0.091504343 0.06203670 -0.20888792  0.035310492
13* y >= 4   1.5  13  59 -0.066840801 0.05926664 -0.17901188  0.051240296
14* y >= 4   1.5  14  59 -0.041564743 0.05922627 -0.15273067  0.077362541
15* y >= 4   1.5  15  59 -0.017079101 0.06238757 -0.13676577  0.104427505
16* y >= 4   1.5  16  59  0.005213193 0.06764387 -0.12999943  0.132730499
17* y >= 4   1.5  17  59  0.023909209 0.07307822 -0.11465525  0.171848255
18* y >= 4   1.5  18  59  0.037606015 0.07679322 -0.12173234  0.179289079
19* y >= 4   1.5  19  59  0.045253060 0.07754647 -0.10315546  0.200192158
20* y >= 4   1.5  20  59  0.047209312 0.07566747 -0.09637584  0.194004388
21* y >= 4   1.5  21  59  0.044186120 0.07238835 -0.09685151  0.182210688
22* y >= 4   1.5  22  59  0.036894834 0.06959524 -0.09268993  0.175323057
23* y >= 4   1.5  23  59  0.026046803 0.06965134 -0.10728992  0.162345053
24* y >= 4   1.5  24  59  0.012353376 0.07459724 -0.13440317  0.157112411
25* y >= 4   1.5  25  59 -0.003474099 0.08506170 -0.17360616  0.159039235
26* y >= 4   1.5  26  59 -0.020724271 0.10019078 -0.23582239  0.158513147
27* y >= 4   1.5  27  59 -0.038685793 0.11853621 -0.27319872  0.191824969
28* y >= 4   1.5  28  59 -0.056765873 0.13878805 -0.34821707  0.195101679
           Pr(Contrast>0)
 1  y >= 4         0.0265
 2  y >= 4         0.0118
 3  y >= 4         0.0037
 4* y >= 4         0.0008
 5* y >= 4         0.0005
 6* y >= 4         0.0020
 7* y >= 4         0.0043
 8* y >= 4         0.0078
 9* y >= 4         0.0127
10* y >= 4         0.0235
11  y >= 4         0.0412
12* y >= 4         0.0685
13* y >= 4         0.1278
14* y >= 4         0.2455
15* y >= 4         0.3900
16* y >= 4         0.5315
17* y >= 4         0.6258
18* y >= 4         0.6845
19* y >= 4         0.7145
20* y >= 4         0.7330
21* y >= 4         0.7300
22* y >= 4         0.7018
23* y >= 4         0.6482
24* y >= 4         0.5662
25* y >= 4         0.4892
26* y >= 4         0.4235
27* y >= 4         0.3738
28* y >= 4         0.3412

Redundant contrasts are denoted by *

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

The time-related treatment effect is similar to what was seen in the earlier frequentist PO model fit.

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.

The initial state (baseline ordinal status) is always conditioned on, so is not marginalized out. Calculations can be done separately for each initial state.
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.
Likewise for later days.
Death is an absorbing state so once death occurs all future status values are also set to death.
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. Daily measurements are too course for the median to work properly. 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). Also we can sum such probabilities over days to estimate mean time in state.
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. A more comprehensive approach was used for ORCHID but here we take a simple approach. Considering only two baseline covariates, find the median initial state and for that state find the median age.

w[, .(table(y0), median=median(as.numeric(as.character(y0))))]

   y0     N median
1:  4  3596      5
2:  5 10873      5
3:  6  4449      5
4:  7  6834      5

w[y0 == '5', median(age)]    # 59

Age 
[1] 59

covar <- list(age=59, yprev='5')

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

The R Hmisc package soprobMarkovOrdm function makes this easy for a single covariate value combination. It computes all the state occupancy probabilities for all posterior draws (here, 4000 of them). From this one can derive any quantity of interest. By computing differences in the 4000 probabilities for each time, state, and treatment, we can compute posterior probabilities of differences in state probabilities. Before moving on to difference in mean number of sick days (dead or in hospital for a reason other than infection control), we get, one placebo, transition probabilities are select days then for all days so that we can see whether time interpolation is trustworthy.

covar$tx <- 'Placebo'
up <- function(data, times) soprobMarkovOrdm(bppodc, data, times, as.character(1:8), absorb='8',
                                             tvarname='day', pvarname='yprev')
u <- up(covar, 1:28)
# Compute posterior mean of probs
p1 <- apply(u, 2:3, mean)
round(p1, 2)

      1    2    3    4    5    6    7    8
1  0.00 0.01 0.00 0.12 0.74 0.12 0.01 0.00
2  0.00 0.03 0.01 0.16 0.59 0.15 0.04 0.00
3  0.00 0.07 0.02 0.18 0.50 0.15 0.08 0.00
4  0.01 0.10 0.02 0.18 0.44 0.14 0.10 0.00
5  0.02 0.14 0.02 0.18 0.39 0.12 0.12 0.01
6  0.03 0.17 0.02 0.18 0.34 0.10 0.13 0.01
7  0.05 0.21 0.02 0.18 0.30 0.09 0.14 0.01
8  0.06 0.24 0.03 0.17 0.27 0.08 0.14 0.02
9  0.08 0.27 0.03 0.16 0.24 0.07 0.14 0.02
10 0.10 0.30 0.03 0.15 0.21 0.06 0.14 0.02
11 0.12 0.32 0.02 0.14 0.19 0.05 0.13 0.02
12 0.14 0.34 0.02 0.13 0.17 0.05 0.13 0.03
13 0.15 0.36 0.02 0.12 0.15 0.05 0.12 0.03
14 0.16 0.38 0.02 0.11 0.14 0.04 0.12 0.03
15 0.17 0.39 0.02 0.10 0.13 0.04 0.11 0.03
16 0.18 0.41 0.02 0.10 0.12 0.04 0.11 0.03
17 0.18 0.42 0.02 0.10 0.11 0.04 0.10 0.04
18 0.18 0.43 0.02 0.09 0.11 0.03 0.10 0.04
19 0.18 0.44 0.02 0.09 0.10 0.03 0.09 0.04
20 0.18 0.45 0.02 0.09 0.10 0.03 0.09 0.04
21 0.18 0.46 0.02 0.09 0.10 0.03 0.09 0.04
22 0.18 0.46 0.02 0.09 0.09 0.03 0.08 0.04
23 0.19 0.47 0.02 0.08 0.09 0.03 0.08 0.05
24 0.19 0.47 0.02 0.08 0.09 0.03 0.07 0.05
25 0.19 0.48 0.02 0.08 0.08 0.02 0.07 0.05
26 0.20 0.48 0.02 0.08 0.08 0.02 0.07 0.05
27 0.21 0.48 0.02 0.08 0.08 0.02 0.06 0.05
28 0.21 0.49 0.02 0.08 0.07 0.02 0.06 0.05

covar$tx <- 'Remdesivir'
u <- up(covar, 1:28)
p2 <- apply(u, 2:3, mean)
round(p2, 2)

      1    2    3    4    5    6    7    8
1  0.00 0.01 0.01 0.13 0.73 0.11 0.01 0.00
2  0.00 0.04 0.02 0.18 0.59 0.13 0.04 0.00
3  0.01 0.08 0.02 0.20 0.49 0.13 0.06 0.00
4  0.02 0.13 0.02 0.21 0.42 0.12 0.08 0.00
5  0.03 0.17 0.03 0.20 0.37 0.10 0.09 0.00
6  0.05 0.22 0.03 0.20 0.32 0.09 0.10 0.01
7  0.07 0.26 0.03 0.19 0.28 0.07 0.10 0.01
8  0.10 0.29 0.03 0.17 0.24 0.06 0.10 0.01
9  0.12 0.32 0.03 0.16 0.20 0.05 0.10 0.01
10 0.15 0.35 0.03 0.14 0.18 0.05 0.10 0.01
11 0.17 0.37 0.02 0.13 0.16 0.04 0.09 0.02
12 0.19 0.39 0.02 0.12 0.14 0.04 0.09 0.02
13 0.20 0.40 0.02 0.11 0.12 0.03 0.09 0.02
14 0.21 0.42 0.02 0.10 0.11 0.03 0.08 0.02
15 0.21 0.44 0.02 0.09 0.11 0.03 0.08 0.02
16 0.21 0.45 0.02 0.09 0.10 0.03 0.08 0.02
17 0.21 0.46 0.02 0.09 0.09 0.03 0.07 0.02
18 0.20 0.48 0.02 0.09 0.09 0.03 0.07 0.03
19 0.20 0.48 0.02 0.08 0.09 0.03 0.07 0.03
20 0.20 0.49 0.02 0.08 0.09 0.02 0.07 0.03
21 0.20 0.50 0.02 0.08 0.09 0.02 0.06 0.03
22 0.20 0.50 0.02 0.08 0.08 0.02 0.06 0.03
23 0.20 0.50 0.02 0.08 0.08 0.02 0.06 0.03
24 0.20 0.50 0.02 0.08 0.08 0.02 0.06 0.03
25 0.21 0.51 0.02 0.08 0.08 0.02 0.05 0.03
26 0.21 0.51 0.02 0.08 0.07 0.02 0.05 0.03
27 0.22 0.51 0.02 0.08 0.07 0.02 0.05 0.03
28 0.23 0.51 0.02 0.08 0.07 0.02 0.04 0.03

cat('\nPlacebo - Remdesivir State Occupancy Probabilities\n\n')


Placebo - Remdesivir State Occupancy Probabilities

round(p1 - p2, 2)

       1     2 3     4    5    6    7    8
1   0.00  0.00 0 -0.02 0.00 0.02 0.00 0.00
2   0.00 -0.01 0 -0.02 0.01 0.02 0.01 0.00
3   0.00 -0.02 0 -0.02 0.01 0.02 0.02 0.00
4  -0.01 -0.03 0 -0.02 0.01 0.02 0.02 0.00
5  -0.01 -0.04 0 -0.02 0.02 0.02 0.03 0.00
6  -0.02 -0.04 0 -0.01 0.02 0.02 0.03 0.00
7  -0.02 -0.05 0 -0.01 0.03 0.02 0.04 0.00
8  -0.03 -0.05 0  0.00 0.03 0.02 0.04 0.01
9  -0.04 -0.05 0  0.00 0.03 0.02 0.04 0.01
10 -0.04 -0.05 0  0.01 0.03 0.01 0.04 0.01
11 -0.05 -0.05 0  0.01 0.03 0.01 0.04 0.01
12 -0.05 -0.05 0  0.01 0.03 0.01 0.04 0.01
13 -0.05 -0.05 0  0.01 0.03 0.01 0.04 0.01
14 -0.05 -0.04 0  0.01 0.02 0.01 0.03 0.01
15 -0.04 -0.05 0  0.01 0.02 0.01 0.03 0.01
16 -0.04 -0.05 0  0.01 0.02 0.01 0.03 0.01
17 -0.03 -0.05 0  0.01 0.02 0.01 0.03 0.01
18 -0.02 -0.05 0  0.01 0.02 0.01 0.03 0.01
19 -0.02 -0.04 0  0.01 0.01 0.01 0.03 0.01
20 -0.02 -0.04 0  0.00 0.01 0.01 0.02 0.01
21 -0.01 -0.04 0  0.00 0.01 0.01 0.02 0.01
22 -0.01 -0.04 0  0.00 0.01 0.01 0.02 0.01
23 -0.01 -0.04 0  0.00 0.01 0.00 0.02 0.01
24 -0.01 -0.03 0  0.00 0.01 0.00 0.02 0.01
25 -0.01 -0.03 0  0.00 0.01 0.00 0.02 0.01
26 -0.01 -0.03 0  0.00 0.01 0.00 0.02 0.02
27 -0.02 -0.02 0  0.00 0.01 0.00 0.01 0.02
28 -0.02 -0.02 0  0.00 0.01 0.00 0.01 0.02

Plot the unconditional probability of being at home and the probability of being dead as a function of time, by treatment.

r <- rbind(data.frame(tx='Placebo',    day=1:28, status='home', p=p1[, 1] + p1[, 2]),
           data.frame(tx='Remdesivir', day=1:28, status='home', p=p2[, 1] + p2[, 2]),
           data.frame(tx='Placebo',    day=1:28, status='dead', p=p1[, 8]),
           data.frame(tx='Remdesivir', day=1:28, status='dead', p=p2[, 8]))
ggplot(r, aes(x=day, y=p, color=tx, linetype=status)) + geom_line() +
  xlab('Day') + ylab('Probability')

Markov Model Without Covariate Adjustment

Let’s repeat the exercise with no covariate adjustment to check how we are estimating the proportion of patients at home at each observed day.

First estimate occupancy probabilities empirically for those in initial state y=5 (the most prevalent one). Death is carried forward.

w <- a[y0 == 5, .(dlast=day, last8 = day == max(day) & y == 8, tx=tx, dlasthosp, dlasthospf), by=id]
w <- w[last8 == TRUE, .(day = (dlast + 1) : 28, y=rep(8, 28 - dlast), tx=tx, dlasthosp, dlasthospf), by=id]
label(w$day) <- label(d$day)
label(w$y)   <- label(d$y)
w <- rbind(a[, .(id, day, y, tx, dlasthosp, dlasthospf)], w)
setkey(w, id, day)
dcf <- w

s <- dcf[, .(phome = mean(y < 3)), by=.(tx, day)]
s[, type := 'Empirical']
empirical <- s
s$tx <- as.character(s$tx)
s$day <- as.integer(s$day)

Now fit the model with which to compare the empirical proportions.

bplain <- blrm(Ocens(yl, yu) ~ yprev + rcs(day, 4) * tx,
               ~ day, cppo = function(y) y,
               data=a, file='bplain.rds', refresh=100)

tlev  <- unique(a$tx)

for(tx in tlev) {   # separately for treatments
    covar$tx <- tx
    # covar still has yprev='5'
  u <- soprobMarkovOrdm(bplain, covar, 1:28, as.character(1:8), absorb='8',
                                             tvarname='day', pvarname='yprev')
  # Compute posterior mean of probs
  p <- apply(u, 2:3, mean)
  s <- rbind(s, data.frame(tx, day=1:28, phome=p[, 1] + p[, 2],
                                                 type='Model'), fill=TRUE)
}
ggplot(s, aes(x=day, y=phome, color=tx, linetype=type)) + geom_line() +
  xlab('Day') + ylab('Probability of Being at Home') +
    labs(color='', linetype='')

Bayesian Inference for Mean Time Not Recovered

We go back to the age-adjusted partial proportional odds Markov model and estimate for each treatment the posterior distribution of mean time not recovered, which is the mean time with y > 3 (which includes death). These are estimated at the median age (59) and modal initial state (y=5).

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

mt <- list()
for(tx in tlev) {
    covar$tx <- tx
  u <- soprobMarkovOrdm(bppodc, covar, 1:28, as.character(1:8), absorb='8',
                                             tvarname='day', pvarname='yprev')
  # Compute 4000 mean number of days with y > 3
  mtime <- function(u) sum(u[, 4 : 8])
  mt[[tx]] <- apply(u, 1, mtime)
}

Compute by treatment the posterior mean expected and 0.95 highest posterior density interval for days not recovered, and for the difference in treatments. Show the posterior distribution for the difference, and compute the probability that drug is better than placebo.

g(mt$Placebo)

    mean    lower    upper 
14.60173 13.19856 16.06856

g(mt$Remdesivir)

    mean    lower    upper 
12.88650 11.48769 14.25641

g(delta <- mt$Placebo - mt$Remdesivir)

     mean     lower     upper 
1.7152342 0.6111744 2.9219856

plot(density(delta), xlab='Placebo - Remdesivir Mean Days Not Recovered', main='')
P <- mean
P(delta > 0)

[1] 0.9985

Let’s repeat the calculations for four initial states and five ages (approximately the quintiles).

donecovs <- FALSE
if(donecovs) {
  w <- readRDS('mt.rds')
  dens <- w$dens
  r    <- w$r
} else {
r <- dens <- NULL
for(yprev in c('4', '5', '6', '7')) {
  for(age in c(45, 55, 65, 75)) {
    cat(yprev, age, '\n')
    mt <- list()
    for(tx in tlev) {
      covar <- data.frame(yprev, age, tx)
      u <- soprobMarkovOrdm(bppodc, covar, 1:28, as.character(1:8), absorb='8',
                            tvarname='day', pvarname='yprev')
      mt[[tx]] <- apply(u, 1, mtime)
    }
    mp <- mt$Placebo
    mr <- mt$Remdesivir
    delta <- mp - mr
    dens <- rbind(dens, data.frame(yprev, age, delta))
    h    <- g(delta)
    r    <- rbind(r,    data.frame(yprev, age, Mean.Placebo=mean(mp), Mean.Remdesivir=mean(mr),
                                   Difference=h[1], Lower=h[2], Upper=h[3], P=P(delta > 0)))
  }
}
rownames(r) <- NULL
w <- list(dens=dens, r=r)
saveRDS(w, 'mt.rds')
}

ggplot(dens, aes(x=delta)) + geom_density() + facet_grid(yprev ~ age) +
  xlab('Difference in Mean Time Not Recovered')
print(r, digits=3)

   yprev age Mean.Placebo Mean.Remdesivir Difference  Lower Upper     P
1      4  45         8.38            7.00      1.381 0.5023  2.37 0.999
2      4  55         9.58            8.06      1.521 0.5407  2.59 0.998
3      4  65        11.76           10.04      1.721 0.5986  2.93 0.998
4      4  75        13.98           12.13      1.851 0.6151  3.14 0.998
5      5  45        12.43           10.84      1.584 0.6138  2.74 0.999
6      5  55        13.75           12.08      1.674 0.5660  2.82 0.999
7      5  65        16.03           14.27      1.760 0.6079  2.99 0.998
8      5  75        18.18           16.43      1.758 0.5639  2.98 0.998
9      6  45        17.72           16.01      1.716 0.5997  2.87 0.998
10     6  55        19.06           17.37      1.689 0.5849  2.85 0.998
11     6  65        21.13           19.56      1.569 0.5619  2.72 0.998
12     6  75        22.85           21.46      1.385 0.4690  2.40 0.998
13     7  45        23.09           22.14      0.945 0.1706  1.71 0.992
14     7  55        23.92           23.07      0.850 0.1198  1.54 0.990
15     7  65        25.09           24.41      0.687 0.1049  1.30 0.987
16     7  75        25.96           25.42      0.538 0.0459  1.01 0.986

Computing Environment

 R version 4.1.0 (2021-05-18)
 Platform: x86_64-pc-linux-gnu (64-bit)
 Running under: Pop!_OS 20.10
 
 Matrix products: default
 BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
 LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
 
 attached base packages:
 [1] splines   stats4    stats     graphics  grDevices utils     datasets 
 [8] methods   base     
 
 other attached packages:
  [1] rstan_2.19.3         StanHeaders_2.21.0-1 VGAM_1.1-3          
  [4] plotly_4.9.3         data.table_1.13.0    rmsb_0.0.2          
  [7] rms_6.2-1            SparseM_1.78         Hmisc_4.6-0         
 [10] ggplot2_3.3.2        Formula_1.2-3        survival_3.2-11     
 [13] lattice_0.20-44

To cite R in publications use:

R Core Team (2021). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.

To cite the Hmisc package in publications use:

Harrell Jr F (2021). Hmisc: Harrell Miscellaneous. R package version 4.6-0, https://hbiostat.org/R/Hmisc/.

To cite the rms package in publications use:

Harrell Jr FE (2021). rms: Regression Modeling Strategies. https://hbiostat.org/R/rms/, https://github.com/harrelfe/rms.

To cite the rmsb package in publications use:

Harrell F (2021). rmsb: Bayesian Regression Modeling Strategies. R package version 0.0.2, https://hbiostat.org/R/rmsb/.

Longitudinal Ordinal Analysis for ACTT-1

Longitudinal Ordinal Analysis for ACTT-1

Introduction

Descriptive Statistics

Data Collection Intensity While in Hospital

Summary of Successive State Transitions

Category Proportions Over Time

Single-Day Analysis

Time to Recovery Analysis with Cox Model

Serial Conditional Discrete Time Markov Ordinal Model

Proportional Odds Assumption

Computing Unconditional Estimates

Markov Model Without Covariate Adjustment

Bayesian Inference for Mean Time Not Recovered

Computing Environment