---
title: Simulating Operating Characteristics of a Sequential Longitudinal Markov Ordinal Outcome Trial
author: |
  | Frank Harrell
  | Department of Biostatistics
  | Vanderbilt University School of Medicine
date: "`r Sys.Date()`"
output:
  rmdformats::readthedown:
    thumbnails: no
    lightbox: no
    gallery: no
    highlight: tango
    use_bookdown: yes
    toc_depth: 3
    fig_caption: yes
    code_folding: show
    embed_fonts: no
    keep_md: false
  pdf_document:
    latex_engine: lualatex
    toc: true
    toc_depth: 3
    fig_caption: yes
    keep_tex: true
    extra_dependencies: ["float"]
urlcolor: blue
linkcolor: blue
header-includes: \usepackage{longtable}
link-citations: yes
description: Markov Simulation
---
<!-- To compile: Rstudio knit to rmdreadthedown or to pdf or knitrg sim2 -->
<style>
.table {
	width: 50%;
}
p.caption {
  font-size: 0.8em;
  color: DarkBlue;
}
</style>
# Summary

The performance of pairwise treatment comparison tests is simulated for 3-level longitudinal ordinal responses.  The setup is a special case of longitudinal continuous response assessments with an override for deaths, which are counted as the highest response category.  Here the continuous daily assessment is dichotomized so that on a given day the response Y=0 for fully recovered and alive, Y=1 for alive but not fully recovered, and Y=2 for dead.  Simulations are done for estimating frequentist test power for both a longitudinal discrete time Markov proportional odds ordinal logistic model and a Cox model for time until recovery, the latter considering death at any time to represent no recovery by 28 days (time to recovery right censored at 28d).  Y=0 and 2 are absorbing states, so once a patient recovers she is assumed to stay recovered.  Under ideal conditions (constant treatment effect over time, proportional odds holds for the treatment effect, proportional hazards holds for time to recovery) the transition model odds ratio equals the Cox model hazard ratio, and the two tests provide identical test statistics and power.  Note that for the general case where there are many levels of Y, the power of the longitudinal ordinal model will greatly exceed that of the time to recovery Cox test as demonstrated [here](https://hbiostat.org/proj/covid19/violet2.html).  Note also that when Y has 3 levels and two of them are absorbing, the Markov transition model has no need to use the previous state as a covariate (since it's constant), which is another reason that the Cox-ordinal longitudinal model equivalence will not hold in general.  For the current limited-Y case the Markov model provides added flexibility and a simpler interpretation in the presence of competing events, by estimating these quantities:

* the probability that Y $\geq 1$ as a function of time and treatment (being not recovered or dead)
* the probability that Y=1 as a function of time and treatment (being alive but not recovered)
* the probability of being dead as a function of time and treatment
* the mean number of days unrecovered (days with Y>0) by treatment

The last quantity is particularly clinically interpretable, being the expected number of days within the 28 day follow-up window for which the patient has died or is alive but not yet recovered.

For comparing two equal size treatment groups of $n=89$ patients each, the frequentist power at $\alpha=0.05$ for both statistical tests is given below, along with the hazard and odds ratios to be detected, and the corresponding reduction in mean days unrecovered.

Hazard Ratio | Reduction in<br>Days Unrecovered | Cox<br>Power
-------------|------------------|-------
1.5          | 2.0              | 0.71
1.6          | 2.3              | 0.82
1.66667      | 2.5              | 0.88
1.75         | 2.7              | 0.93

As shown [below](#ttrdeath), the Cox model analysis for time to recovery in essence makes the proportional odds assumption for the outcomes of recovered, not recovered, and dead.  When the proportional odds assumption is violated, the Cox model hazard ratio for time to recovery is almost exactly equal to the "average" odds ratio from the PO model that falsely assumes proportional odds.  Like the PO model in this case, the Cox model (and its special case the logrank test) is not able to recover a pure effect of a treatment on recovery unless the treatment affects death the same way as it affects recovery (i.e., proportional odds holds).  The same examples show that when the partial PO model is used, by design it is able to accurately recover the separate recovery and mortality treatment effects.


# Overview

Consider a clinical trial for the treatment of Stevens-Johnson syndrome / toxic epidermal necrolysis, a rare and life-threatening severe cutaneous adverse reaction induced primarily by drugs.  Our goal is to estimate Bayesian operating characteristics for a two treatment arm comparison with five pre-planned data looks with optional stopping for efficacy, inferiority, or futility.  A first-order discrete time Markov Bayesian partial proportional odds (PO) model is used, with daily patient assessment on days 1, 2, ..., 28 after randomization.  The underlying daily measurement is the percent body surface area for which skin has not re-attached (re-epitheliazation), with death considered as worse than any level of skin detachment.  The most powerful analysis will treat % BSA involvement as a continuous ordinal measurement with death override.  But we now concern ourselves with a simpler daily outcome: Y=0 (alive and completely recovered), Y=1 (alive but some skin remains detached), and Y=2 (dead).  Consider Y=0 and Y=2 to be absorbing states, so per-patient serial data records stop once either of these is achieved.

Once Markov transition probabilities are transformed into state occupancy probabilities (SOPs), these can be used as primary estimands or a combination of the 28 days of SOPs can be summarized.  In the study the primary estimand is the difference between treatments in the mean time not recovered, which is the difference in the sums over the 28 days of the probability that Y>0.  Were there to be no deaths, the mean time not recovered is just one day less than the mean time to recovery.  For the purpose of the simulations presented here the estimand will be the transition odds ratio in the Markov PO model.  Note that the PO model places death as worse than non-recovery, but makes no assumption about the spacing of categories, i.e., going from death to not recovered is not assumed to be the same as going from not recovered to recovered.

It is common practice to analyze such data by collapsing each patient's raw data into a single number---time to recovery---and then applying a Cox model comparison or its competing risk analog.   For these analyses, death is considered a censoring event, and patients who die are censored by treating them as non-recovered on day 28.  Thus a death at day 2 is considered as serious as a death at day 27 and these are considered as serious as being alive without recovery on day 28.  It can be difficult to interpret time until a good outcome when the outcome can be interrupted by a bad outcome.  State transition models provide estimates of probabilities of easy-to-interpret conditions, e.g. day-specific P(dead), P(alive but not recovered), or P(either being dead or being alive but not recovered).  The easy generalization of such models to fully use the continuous version of the daily outcome (% body surface area involvement) is another advantage.

As an aside, it is customary to summarize death-censored time-to-recovery using the median.  The median requires one to have continuous measurements, and ties in the day of recovery across patients makes the median simultaneously inefficient and non-robust.  Adding one patient to the study can alter the median by a whole day.  One would need to have hourly assessments for the median number of days to have good statistical behavior.  These same problems do not exist for the mean time to recovery.

More statistical background and statistical computing details are available [here](https://hbiostat.org/R/Hmisc/markov/sim.html).  As detailed [here](https://hbiostat.org/proj/covid19/violet2.html), the compound symmetric correlation pattern assumption of a simpler random effects ordinal model is severely violated for designs like ours, hence the use of the Markov serial correlation pattern.  Note that even when the proportional odds model is strongly violated for treatment, a simple function of the odds ratio (OR) from the model almost exactly [reproduces the Wilcoxon statistic](https://www.fharrell.com/post/po) or *concordance probability*.  The Wilcoxon statistic, when scaled to 0-1 (i.e., the concordance probability) is to within an average absolute error of 0.002 equal to $\frac{\mathrm{OR}^{0.66}}{1 + \mathrm{OR}^{0.66}}$.  The most important non-PO effect to model is due to the possible change in the mix of outcomes over follow-up time---an effect not necessarily involving treatment.  This is modeled using the [Peterson and Harrell (1990)](https://rss.onlinelibrary.wiley.com/doi/abs/10.2307/2347760) *partial proportional odds model*.

```{r setlooks}
olooks   <- c(30, 54, 78, 102, 126, 150, 174, 198, 222, 246, 267)
kl       <- length(olooks)
colooks1 <- paste(olooks[- kl], collapse=', ')
colooks  <- paste0(colooks1, ', and ', olooks[kl])
looks    <- olooks * 2 / 3
clooks1  <- paste(looks[- kl], collapse=', ')
clooks   <- paste0(clooks1, ', and ', looks[kl])
```
  
The new study in question is a 3-arm study with a planned maximum sample size of 267 patients with looks after follow-up is complete on `r colooks` (the target overall sample size) patients.  Looks are scheduled after the first 30 patients have completed follow-up, then every 6m over an additional 5y period (23.7 patients per 6m rounded to 24, with 21 enrolled in the final 6m).  Here we consider 2-arm comparisons so the combined sample sizes are `r clooks` (final look at planned maximum sample size).

# R Setup

```{r setup}
outfmt <- if(knitr::is_html_output ()) 'html'  else 'pdf'
markup <- if(knitr::is_latex_output()) 'latex' else 'html'
require(rms)          # automatically engages rms(Hmisc)
require(data.table)   # used to compute various summary measures
if(outfmt == 'html') require(plotly)
fdev <- switch(outfmt, html='png', pdf='pdf')
knitrSet(lang='markdown', w=7, h=7, dev=fdev, fig.path=paste0(fdev, '/sim2-'))
options(prType=markup)
## If producing html, ggplot2 graphics are converted to plotly graphics
## so that hovertext will show extra information
## Assumes aes(..., label=txt) used
ggp <- if(outfmt != 'html')
  function(ggobject, ...) ggobject else Hmisc::ggplotlyr

# When outputting to pdf, make LaTeX place figures right at the point		
# in which they are created; requires float.sty (see yaml header)
if(outfmt == 'pdf') knitr::opts_chunk$set(fig.pos = 'H', out.extra='')

# Function in Hmisc that will sense chunk has results='asis' and
# print in Rmarkdown verbatim mode if needed
pr <- markupSpecs$markdown$pr

# Define the squote, equote, pr, findstart, sop12, and dosim functions
source('funs.r')

# Read the source code for the hashCheck and runifChanged functions from
# https://github.com/harrelfe/rscripts/blob/master/hashCheck.r
# This makes it easy to see if any objects changed that require re-running
# a simulation, and reports on any changes
getRs('hashCheck.r', put='source')

# Load runParallel function from github
# Makes parallel package easier to use and does recombination
# over batches
getRs('runParallel.r', put='source')

# Quick way to get the variance-covariance matrix from a vgam fit, per author Thomas Yee
# See ~/r/VGAM/notes.md
vgamcov <- function(fit) {
  covun <- chol2inv(fit@R)
  cnames <- names(fit@coefficients)
  dimnames(covun) <- list(cnames, cnames)
  covun
}
```

# Simulation Model

To be able to simulate serial 3-level ordinal outcomes (Y=0,1,2 above) it is necessary to have estimates of control arm daily outcome probabilities.  Our estimates come from the corticosteriod arm of the trial by [Wang et al](https://www.jci.org/articles/view/93349), Figure 3 Panel D.  The cumulative incidence estimates of complete re-epitheliazation there correspond to survivors.  The new trial's follow-up period for the primary Y=0,1,2 endpoint is 28 days.  We assume that the 21d mortality for the expected severely ill patient population being studied to be 0.2, and assume an exponential distribution for the time until death such that $S(21) = 0.8$, leading to the cumulative incidence of death by day $t$ of $1 - \exp(\lambda t)$ where $\lambda = \frac{-\log(0.8)}{21} = 0.010625$.  We smooth the recovery probabilities of Wang _et al_ and interpolate the estimates to 1, 2, ..., 28 days, then scale each day's estimate by the probability of the event upon which they are conditioned (being alive).

Wang's estimate for 28d recovery probability (0.615) is too low.  Apply a multiplicative correction to yield 0.8.


```{r controlmodel}
d <- read.csv('wangFig3D.csv', col.names=c('t', 'p'), header=FALSE)
d$p <- d$p / 100
s <- with(d, supsmu(t, p))
w <- approx(s, xout=1:28)
t <- w$x
p <- w$y
lambda <- -log(0.8) / 21
pdeath <- 1 - exp(- lambda * (1 : 28))
precoveredAlive  <- p * (1 - pdeath) * 0.7 / 0.622892705
precoveredAlive2 <- p * (1 - pdeath) * 0.8 / 0.622892705

pnotrecoverdAlive  <- 1 - pdeath - precoveredAlive
pnotrecoverdAlive2 <- 1 - pdeath - precoveredAlive2

type <- rep(c('death', 'alive and not recovered',
		 'alive and recovered with multiplier',
     'dead or not recovered', 'Wang recovery if survive'), each=28)
y    <- c(pdeath, pnotrecoverdAlive2, precoveredAlive2, 1 - precoveredAlive2, p)
g <- ggplot(data.frame(Day=1:28, Probability=y),
            aes(x=Day, y=Probability, color=type)) + geom_line()
ggp(g)
```

Next compute the mean time to recovery in the control arm, and under various hazard ratios.  Power of the ordinary Cox two-sample test is estimated two ways.  First, exponential distributions are assumed so that all that matters is the number of recovered patients.  Second, the Cox test is simulated by drawing recovery time samples from the Wang control arm adjusted to 0.8 final recovery probability, and assuming proportional hazards to get the corresponding active treatment arm simulated recovery times.  Times are simulated by computing the step function inverse of the Wang cumulative distribution function.

The adjusted Wang cumulative recovery curve and the curve that would result from a treatment benefit of HR=1.667 are found below.

```{r curves,bot=1}
s <- 1 - precoveredAlive2
plot(1:28, 1 - s, xlab='', ylab='Cumulative Probability of Recovery', type='l', ylim=c(0, 1))
lines(1:28, y <- 1 - s ^ 1.6666666, col='blue')
title(sub='Black: adjusted estimates from Wang for control\nBlue: estimates for treated group, using HR=1.667', adj=0)
title(sub=paste(round(max(1 - s), 2), 'at 28d\n', round(max(y), 2), 'at 28d'), adj=1)
```

```{r controlhr}
# Approximate power assuming exponential distributions (just need expected # events by tx)
cpow <- function(hr, nevents1, nevents2, alpha=0.05) {
  zcrit <- - qnorm(alpha / 2)
  sd    <- sqrt((1. / nevents1) + (1. / nevents2))
  z     <- abs(log(hr)) / sd
  1. - (pnorm(zcrit - z) - pnorm(- zcrit - z))
}

# Simulated power based on Wang
pval <- function(hr, n1, pl=FALSE) {
  u  <- runif(n1)
  y1 <- approx(c(precoveredAlive2, 1), 1:29, xout=u, method='constant', f=1)$y
  e1 <- ifelse(y1 == 29, 0, 1)
  y1 <- pmin(y1, 28)
  u  <- runif(n1)
  y2 <- approx(c(1. - (1. - precoveredAlive2) ^ hr, 1), 1:29, xout=u, method='constant', f=1)$y
  e2 <- ifelse(y2 == 29, 0, 1)
  y2 <- pmin(y2, 28)
  y  <- c(y1, y2)
  e  <- c(e1, e2)
  tx <- c(rep(0, n1), rep(1, n1))
  if(pl) {
    plot(survfit(Surv(y, e) ~ tx))
    points(1:28, 1. - precoveredAlive2)
  }
  f  <- coxph(Surv(y, e) ~ tx)
  lr <- diff(f$loglik)
  1. - pchisq(2. * lr, 1)
}

# Test the function on large random samples
pval(1.5, 20000, pl=TRUE)

# Also estimate power using Jon Fintzi's delayed exponential model

pvaljf <- function(hr, n1, pl=FALSE) {
  tx <- rep(0:1, each=n1)
  y  <- c(7 + ceiling(rexp(n1, log(5) / (28 - 7))),
          7 + ceiling(rexp(n1, log(5) / (28 - 7) * hr)))
  e <- ifelse(y > 28, 0, 1)
  y[e == 0] <- 28   # non-recovery
  if(pl) {
    plot(survfit(Surv(y, e) ~ tx))
    points(1:28, 1. - precoveredAlive2)
  }
  f <- coxph(Surv(y, e) ~ tx)
  lr <- diff(f$loglik)
  1. - pchisq(2. * lr, 1)
}
pvaljf(1.5, 20000, pl=TRUE)

simpow <- function(hr, n1, fun=pval, nsim=1000) {
  pvals <- numeric(nsim)
  for(i in 1 : nsim) pvals[i] <- fun(hr, n1)
  mean(pvals < 0.05)
}

n1 <- 178 / 2

# Estimate power with JF's methods for HR=1.5
simpow(1.5, n1, pvaljf)

hrs <- c(1, 1.25, 1.5, 1.6, 1.6666667, 1.75, 2)

coxsim <- function() {
s <- 1 - precoveredAlive2
mtr <- pow <- spow <- numeric(0)
for(h in hrs) {
  s2   <- s ^ h
  mtr  <- c(mtr, sum(s2))
  pow  <- c(pow, cpow(h, n1 * max(1. - s), n1 * max(1. - s2)))
  spow <- c(spow, simpow(h, n1, nsim=2000)) 
}
diff <- mtr[1] - mtr
w <- data.frame(HR=hrs, `Mean Time<br>To Recovery`=mtr, Difference=diff,
                Power=pow, `Simulated<br>Power`=spow, check.names=FALSE)
w
}
seed <- 3
set.seed(seed)
w   <- runifChanged(coxsim, seed, n1, hrs, precoveredAlive2, simpow, pval)

knitr::kable(w, digits=3,
   caption='Hazard ratio, mean time to recovery, and power of Cox test')
```

**NOTE**: Too many simulations failed to have the VGAM::vglm fit converge to trust the results with a target cumulative incidence of 0.8.  For what follows 0.7 is used.

Compute parameters of a Markov ordinal model that provide the best compromise fit the the SOPs for death and for alive not recovered, over 28d.   Since there are only 3 states and two of them are absorbing, the previous state is always Y=1 so previous state is ignored in the transition functions `g1` and `g2`.  `g1` and `g2` are the same for control patients.  `g1` parameterizes a constant treatment effect over time, and `g2` parameterizes a linearly increasing effect such that there is no effect at day 1 and the maximum effect is reached at day 10.  The `parameter` for `g1` is the log odds transition ratio for treatment, and for `g2` is the log odds ratio for treatment at 10d.

```{r estmod}
# Define function that computes logit of transition probabilities
# Allow for non-PO with regard to a quadratic time effect

g1 <- function(yprev, t, gap, X, parameter=0, extra) {
  if(! all(yprev == 1)) {prn(yprev); stop('program logic error')}
	kappa <- extra[1:4]
	lp <- matrix(0.,
	             nrow=length(yprev), ncol=2,
	             dimnames=list(as.character(yprev), c('1', '2')))
	# 2 columns = no. distinct y less 1 = length of intercepts
	# lp[yp, ] is a 3-vector because the kappa components are split out into a 3-vector
	# kappa[1:2] are the PO linear and quadratic time trend parameters and
	# kappa[3:4] are the non-PO extra values for death
	for(yp in yprev)
 	  lp[as.character(yp), ] <- kappa[1] * t + kappa[2] * t ^ 2 +
	                            cbind(0, 1) * (kappa[3] * t  + kappa[4] * t ^ 2) +
	                            parameter * (X == 2)
  lp
}

g2 <- function(yprev, t, gap, X, parameter=0, extra) {
  if(! all(yprev == 1)) {prn(yprev); stop('program logic error')}
	kappa <- extra[1:4]
	lp <- matrix(0.,
	             nrow=length(yprev), ncol=2,
	             dimnames=list(as.character(yprev), c('1', '2')))
	# 2 columns = no. distinct y less 1 = length of intercepts
	# lp[yp, ] is a 3-vector because the kappa components are split out into a 3-vector
	# kappa[1:2] are the PO linear and quadratic time trend parameters and
	# kappa[3:4] are the non-PO extra values for death
	for(yp in yprev)
 	  lp[as.character(yp), ] <- kappa[1] * t + kappa[2] * t ^ 2 +
	                            cbind(0, 1) * (kappa[3] * t  + kappa[4] * t ^ 2) +
	                            parameter * (X == 2) * pmin(t - 1, 9) / 9
  lp
}


# SOP targets for control patients:
target <- cbind(precoveredAlive, pnotrecoverdAlive, pdeath)
rownames(target) <- as.character(1:28)

seed      <- 6
int.try   <- c(12, -4.4)
extra.try <- c(-1, .02, .9, -.02)
fun       <- function(...) findstart(intercepts=int.try, extra=extra.try,
  	  					             initial=1, absorb=c(0, 2),
	  	  				             times=1:28, y=0:2, g=g1, seed=seed)

sparams   <- runifChanged(fun, seed, target, g1, int.try, extra.try, findstart, intMarkovOrd,
                          file='sparams.rds', .inclfun.=FALSE)

sop12(y=0:2, times=1:28, initial=1, absorb=c(0,2),
      intercepts=sparams$intercepts, g=g1, extra=sparams$extra, which=1)

# Save the optimal values of the extra vector as default values for the extra argument
# so we don't need to specify them in later steps
formals(g1)$extra <- sparams$extra
formals(g2)$extra <- sparams$extra
saveRDS(list(g1=g1, g2=g2), 'gsim2.rds')

# Save intercepts
ints <- sparams$intercepts
```

Parameter values were found that made the Markov PO model accurately reflect the existing data, with a total absolute error of 0.59 over 84 state occupancy probabilities.  Note that a strong non-proportional odds component for the absolute follow-up time effect was necessary to reflect the changing relative frequencies of the three outcomes over time.

# Transition Odds Ratios vs. Change in Mean Time Unrecovered

The primary treatment effect in a Markov PO model is stated in terms of the transition odds ratio (OR), i.e., the ratio of odds of $Y(t) \geq y | Y(t-1) = a$ for treatment B vs. treatment A (with only one non-absorbing state the $Y(t-1) = a$ condition is superfluous).  The primary result will be stated in terms of the B - A difference in mean time unhealed (which includes time not alive).  We need to be able to relate one to the other.  Two simulations are done.  In the first, the treatment effect is constant over time in the simulated data and in the model used to analyze them.  In the second simulation we allow in the simulations and in the final analysis for a treatment $\times$ time interaction.  For these simulations we define the OR $\exp{\beta}$ to be the treatment effect at 10d, assume the OR stays there after 10d, and assume that the log odds ratio is linear in time from $t=1$ to $t=10$ with the OR being 1.0 at $t=1$ (no treatment effect on the first day post randomization).  So define the treatment effect on day $t$ on the log odds scale to be $\frac{\min(t - 1, 9)}{9} \beta$.  But for the fitted model for these second simulations we use a quadratic time $\times$ treatment interaction to reflect imperfect knowledge about a time-varying treatment effect.

In addition to computation of days unrecovered from the true model, simulate a sample of 20,000 patients for each model and report the mean and median time to recovery.  Note that because of censoring, the median is not valid when its estimate exceeds 27d.   It is also very problematic because the variable being summarized is not continuous.  The mean does not have a problem with discrete data.

Vary treatment effect $\beta$ and compute the difference and ratio of mean days not recovered, separately for the two true treatment effect models.  Note that time to recovery is defined so that it is expected to be approximately one day longer than time unrecovered.

```{r orvsmt,results='asis'}
seed <- 17
set.seed(seed)
ors <- sort(c(1/1.6, seq(0.1, 1, by=0.05)))
m <- length(ors)
# Define time to recovery with treatment of death at any time as censored at 28d
tcl <- function(y) {
  if(any(y == 2)) return(list(rtime=28L, event=0L))   # death
  if(any(y == 0)) return(list(rtime=min(which(y == 0)), event=1L))
	# which works since time=1,2,3,...
  list(rtime=length(y), event=0L)     # censored; length(y) should be 28
}
fun <- function() {
  W <- list()
  for(wg in 1 : 2) {
    mnr <- ttrmean <- ttrmedian <- numeric(m)
    w <- NULL
    i   <- 0
    for(or in ors) {
      gg <- switch(wg, g1, g2)
      sop <- soprobMarkovOrd(0:2, 1:28, initial=1, absorb=c(0,2),
		                         intercepts=ints, g=gg,
                             X=2, parameter=log(or))
    i <- i + 1
    mnr[i] <- sum(1. - sop[, 1])   # mean time not recovered
    s <- simMarkovOrd(n=20000, y=0:2, 1:28, initial=1, X=c(group=2),
	                    absorb=c(0, 2), intercepts=ints, g=gg, parameter=log(or))
	  setDT(s, key='id')
	  ts <- s[, tcl(y), by=id]
	  ttrmean  [i] <- mean(ts$rtime)
	  ttrmedian[i] <- median(ts$rtime)
    }
    W[[wg]] <- 
      data.frame(OR=ors,
                 `Mean Days Unrecovered` = mnr,
                 `Difference in Days`    = mnr[m] - mnr,
	  			  			Ratio                  = mnr / mnr[m], 
	  				  		'Mean TTR'             = ttrmean,
	  					  	'Median TTR'           = ttrmedian,
	  						  check.names=FALSE)
  }
  W
}

W <- runifChanged(fun, seed, ors, tcl, g1, g2, ints, .inclfun.=FALSE)

for(wg in 1 : 2) {
  cap <- switch(wg, 'Constant treatment effect', 'Changing treatment effect')
  print(knitr::kable(W[[wg]], digits=3, caption=cap))
  }
```

We may focus on a transition OR=0.6 which translates to a 2.6 reduction in days not recovered (2.4d for time-varying treatment effect model) or ratio of mean times of 0.88 and 0.80 or a hazard ratio of 1/0.6 = 1.667.

# Effect of OR=0.6 on Time to Recovery Distribution

From our model based on Wang et al's result (after multiplicative factor applied), derive the time to recovery distribution under the null and under OR=0.6 by drawing two samples of size n=20,000.

```{r hreffect}
set.seed(5)
s1 <- simMarkovOrd(n=20000, y=0:2, 1:28, initial=1, X=c(group=1),
                  absorb=c(0, 2), intercepts=ints, g=g1, parameter=0.)
s2 <- simMarkovOrd(n=20000, y=0:2, 1:28, initial=1, X=c(group=2),
	                 absorb=c(0, 2), intercepts=ints, g=g1, parameter=log(0.6))
s  <- rbind(s1, s2)
s$group <- factor(s$group, 1:2, c('control', 'active'))
setDT(s, key=c('group','id'))
ts <- s[, tcl(y), by=.(group, id)]
f  <- npsurv(Surv(rtime, event) ~ group, data=ts)
w  <- round(1. - summary(f, times=28)$surv, 3)
survplot(f, fun=function(y) 1. - y, ylab='Cumulative Probability of Recovery')
```

A OR0.6 moves a cumulative probability of recovery by 28d from `r w[1]` to `r w[2]`.

<a name="ttrdeath"></a>

# Effect of Death on Time to Recovery

The interpretation of Cox model parameters under censoring of time to recovery at 28d for death at any time is worth exploring.  Under the proportional odds (PO) assumption, a treatment that affects mortality to the same degree that it affects non-recovered will have fewer deaths and a shorter time to recovery.  If there is non-PO to the extent that there is no treatment effect on mortality in the Markov PO model, deaths are random events, and even though the time to recovery is censored at 28d, the random deaths may be expected to weaken the treatment effect.  Put that to the test by simulating for each degree of non-PO a 20000 patient dataset with 10000 in each treatment group.  Estimate the hazard ratio for recovery for each case and compare its reciprocal to to the group 2:1 OR for $Y\geq 1$ and the OR for $Y=2$.  The main OR (for $Y \geq 1$) is 0.75 in all cases.  The OR estimated from a PO Markov model and from a partial PO model are also presented.  PPO in the table refers to a fitted partial PO model that unlinks the death effect from the recovery effect.

```{r nonpottr}
file <- 'orhrdeath.rds'
require(VGAM)
seed <- 11
set.seed(seed)
n <- 20000
g3 <- function(yprev, t, gap, X, parameter=0, extra) {
  if(! all(yprev == 1)) stop('program logic error')
	kappa <- extra[1:4]
	tau   <- extra[5]   # amount of non-PO for treatment for Y=2
	lp <- matrix(0.,
	             nrow=length(yprev), ncol=2,
	             dimnames=list(as.character(yprev), c('1', '2')))
	# 2 columns = no. distinct y less 1 = length of intercepts
	# lp[yp, ] is a 3-vector because the kappa components are split out into a 3-vector
	# kappa[1:2] are the PO linear and quadratic time trend parameters and
	# kappa[3:4] are the non-PO extra values for death
	for(yp in yprev)
 	  lp[as.character(yp), ] <- kappa[1] * t + kappa[2] * t ^ 2 +
	                            cbind(0, 1) * (kappa[3] * t  + kappa[4] * t ^ 2) +
	                            (X == 2) * cbind(parameter, parameter + tau)
  lp
}

formals(g3)$extra <- c(sparams$extra, 0.)

hashobj <- hashCheck(seed, n, g3, ints, file=file)
hash    <- hashobj$hash
z       <- hashobj$result
if(! length(z)) {
  s1 <- simMarkovOrd(n=n / 2, y=0:2, 1:28, initial=1, X=c(group=1),
     								 absorb=c(0, 2), intercepts=ints, g=g3, parameter=0)
  setDT(s1, key='id')
  t1 <- s1[, tcl(y), by=id]

  for(tau in c(-0.3, 0, 0.3, 0.6, 0.9, 1.5)) {
    formals(g3)$extra <- c(sparams$extra, tau)
    # tau ignored when X=1
    s2 <- simMarkovOrd(n=n / 2, y=0:2, 1:28, initial=1, X=c(group=2),
        							 absorb=c(0, 2), intercepts=ints, g=g3, parameter=log(0.75))
    setDT(s2, key='id')
    s12 <- rbind(s1, s2)
    s12[, group := factor(group)]
    f <- try(vglm(y ~ pol(time, 2) + group, cumulative(parallel=FALSE ~ pol(time, 2), reverse=TRUE),
              data=s12), silent=TRUE)
    or.overall <- if(inherits(f, 'try-error')) NA else exp(coef(f)['group2'])
    f <- try(vglm(y ~ pol(time, 2) + group, cumulative(parallel=FALSE ~ pol(time, 2) + group, reverse=TRUE),
              data=s12), silent=TRUE)
    if(inherits(f, 'try-error')) or1 <- or2 <- NA else {
      or1 <- exp(coef(f)['group2:1'])
      or2 <- exp(coef(f)['group2:2'])
    }
  
	  t2 <- s2[, tcl(y), by=id]
	  w  <- rbind(cbind(t1, group='1'), cbind(t2, group='2'))
	  f  <- survival::coxph(Surv(rtime, event) ~ group, data=w)
	  z  <- rbind(z, data.frame(tau=tau,
	                            'True OR for<br>not recovered'=0.75,
	                            'True OR<br>for death'=exp(log(0.75) + tau),
	                            'Estimated OR<br>under PO'=or.overall,
	                            'PPO Estimated OR<br>for not recovered'=or1,
	                            'PPO Estimated OR<br>for death'=or2,
	                            '1/HR'=round(exp(-coef(f)), 3),
	                            check.names=FALSE))
  }
  setattr(z, 'hash', hash)
  saveRDS(z, file=file, compress='xz')
}
knitr::kable(z, digits=3, row.names=FALSE)
```

See that when the treatment effect on death is in the wrong direction (tau > -log(0.75)), the reciprocal of the Cox model estimate of the hazard ratio for time to recovery is moved towards the null and away from the effect of treatment on recovery per se.  When the treatment benefits mortality more than it benefits recovery, the Cox hazard ratio is more impressive than the OR for recovery.  So analysis of time to recovery is mixing in the mortality effect, similar to the way the PO model does.  One could go almost as far as saying that the competing risk Cox model analysis for time to recovery also assumes PO.  The time to recovery HR (actually its reciprocal, to convert to a "larger is bad" framework) is very close to the average OR for treatment that falsely assumes PO.

# Clinical Trial Simulation

In addition to simulating longitudinal responses under a Markov partial proportional odds model, the `estSeqMarkovOrd` at the last data look also simulates a Cox proportional hazards $\chi^2$ statistic for time to recovery.  Time to recovery is the first time at which Y=0, and if none or any Y=2, time to recovery is censored at 28d.  The purpose of this additional simulation is to compare the frequentist power of a two-sample comparison of time until recovery with the frequentist power of the transition odds ratio at days 10-28 in the Markov proportional odds model.

## Simulate a Single Large Trial

The simulation model is given above and is encapsulated in the `g` functions previously defined.  For `g2` the analysis model is not informed of the linearly increasing treatment effect from day 1 to day 10 that stays constant from days 10-28.  Instead, the analysis model used on each simulated dataset has a quadratic time $\times$ treatment interaction.

Let's simulate a single 10000-patient trial to check that we can recover a true log odds ratio of log(0.5) = -0.69 overall (model 1) or at $t=10$ (model 2).

```{r simlarge}
require(VGAM)
n <- 10000
for(wg in 1 : 2) {
  cat(if(wg == 1) 'Constant treatment effect\n\n' else '\nTime-varying treatment effect\n\n')
  set.seed(8)
  g <- if(wg == 1) g1 else g2
  s1 <- simMarkovOrd(n=n / 2, y=0:2, 1:28, initial=1, X=c(group=1),
  									absorb=c(0, 2), intercepts=ints, g=g, parameter=0)
  s2 <- simMarkovOrd(n=n / 2, y=0:2, 1:28, initial=1, X=c(group=2),
	  								absorb=c(0, 2), intercepts=ints, g=g, parameter=log(0.5))
  s2$id <- s2$id + n
  s <- rbind(s1, s2)
  s$yprev <- as.factor(s$yprev)
  s$group <- as.factor(s$group)
  form <-  y ~ pol(time, 2) + group
  if(wg == 2) form <- y ~ pol(time, 2) * group
  f <- vglm(form, cumulative(parallel=FALSE ~ pol(time, 2), reverse=TRUE), data=s)
  k <- coef(f)
  # Compute treatment effect; is a constant for model 1
  if(wg == 1) cat('beta =', round(k['group2'], 3), '\n')
  else {
    # Plot estimated treatment effect over time
    time <- 1 : 28
    effect <- k[7] + time * k[8] + time * time * k[9]
    plot(time, effect, type='l')
    lines(c(1, 10, 28), c(0, log(0.5), log(0.5)), col=gray(0.8))
    # Repeat, doing this another way
    contrvgam <- function(fit) {  #, n1=NULL, n2=NULL, parameter=NULL) {   # for customized vgam contrasts if non-PO in effect
	    beta   <- coef(fit)
	    i      <- c('group2', 'pol(time, 2)time:group2', 'pol(time, 2)time^2:group2')
	    b      <- beta[i]
	    effect <- b[1] + 10 * b[2] + 100 * b[3] 
	    # Treatment effect at 10 days (quadratic time x treatment interaction)
	    # if(length(n1)) return(list(Contrast=effect, SE=approxse(or=exp(parameter), n1=n1, n2=n2)))
	    v      <- vgamcov(fit)
	    if(any(i %nin% rownames(v))) {saveRDS(fit, 'problemfit.rds'); stop()}
	    v      <- v[i, i]
	    vr     <-       v[1, 1] +  100 * v[2, 2] + 10000 * v[3, 3] + 20 * v[1, 2] +
	              200 * v[1, 3] + 2000 * v[2, 3]
	    list(Contrast=effect, SE=sqrt(vr))
    }
  print(contrvgam(f))
  }
}
```

The partial PO Markov model was able to recover the treatment effect large simulated data, as the fitted quadratic effect was not far from the true linear spline time-related treatment effect shown in gray above.

## Simulation to Approximate the Standard Error of the Treatment Effect

Let's see if we can satisfactorily approximate the standard error of the treatment effect.  This would speed up simulations in the future.  We model on the scale of the log of the $\hat{\beta}$ standard error.  Do this separately for model 1 (one treatment parameter) and model 2 (three treatment parameters).

```{r approxv}
# Try seeing how predictable is the standard error of the main contrast, for varying true effects and sample sizes
file <- 'vn.rds'
seed <- 1
set.seed(seed)
ors     <- seq(0.1, 1, by=0.1)
combos  <- expand.grid(or=ors, n=c(50, 100, 150, 200, seq(300, 2000, by=100)), model=1:2)
hashobj <- hashCheck(seed, combos, g1, g2, ints, file=file)
hash    <- hashobj$hash
w       <- hashobj$result
if(! length(w)) {
  w <- combos
  v <- numeric(nrow(w))
  for(i in 1 : nrow(w)) {
    n <- w[i, 'n']
    parameter <- log(w[i, 'or'])
    m <- w[i, 'model']
    g <- if(m == 1) g1 else g2
    s1 <- simMarkovOrd(n=n / 2, y=0:2, 1:28, initial=1, X=c(group=1),
	   								absorb=c(0, 2), intercepts=ints, g=g, parameter=0)
    s2 <- simMarkovOrd(n=n / 2, y=0:2, 1:28, initial=1, X=c(group=2),
		  							absorb=c(0, 2), intercepts=ints, g=g, parameter=parameter)
    s <- rbind(s1, s2)
    s$group <- as.factor(s$group)
    form <-  y ~ pol(time, 2) + group
    if(m == 2) form <- y ~ pol(time, 2) * group
    f <- try(vglm(form, cumulative(parallel=FALSE ~ pol(time, 2), reverse=TRUE), data=s), silent=TRUE)
    v[i] <- if(inherits(f, 'try-error')) NA else if(m == 1) vcov(f)['group2', 'group2'] else contrvgam(f)$SE ^ 2
  }
  w$v <- v
  setattr(w, 'hash', hash)
  saveRDS(w, file, compress='xz')
}
ratio <- subset(w, model == 2)$v / subset(w, model == 1)$v
pr('Median model 2 : model 1 ratio of variances of treatment effect', inline=median(ratio))
w <- subset(w, v < 1)   # remove outliers from vgam problems
v <- w$v
ggplot(w, aes(x=n, y=v * n, shape=factor(model), color=or)) + geom_point()
```

Judging by the median variance ratio and the graph above, the 3-parameter treatment contrast is less than half as efficient as the 1-parameter contrast assuming constant treatment effect.  Judging by the flatness of the above trends in the variance of $\hat{\beta}$ multiplied by the sample size, the dominant form of the variance is $\frac{k}{n}$.

```{r fixap}
# Fit various models to log(v)
h <- function(form) {
  z <- NULL
  for(type in c('ols', 'robust')) {
    h <- switch(type,
                robust = MASS::rlm(form, data=dat),
                ols    = ols(form, data=dat))
    r <- resid(h)
    ssr <- var(fitted(h))
    r2 <- ssr / (ssr + var(r))
    z <- rbind(z, data.frame(type=type, model=as.character(form)[3], 
                             stat=c('Mean', 'Median', 'P(|error| > 0.1)'),
                             error=c(mean(abs(r)), Median=median(abs(r)),
                               mean(abs(r) > 0.1)),
                             r2=r2,
                             check.names=FALSE))
  }
  z
}

for(wg in 1 : 2) {
  cat('\nModel', wg, '\n\n')
  dat <- subset(w, model == wg)
  z <- rbind(
    h(log(v) ~ log(n)),
    h(log(v) ~ pol(log(n), 2)),
    h(log(v) ~ log(n) + log(or)),
    h(log(v) ~ log(n) * log(or)),
    h(log(v) ~ pol(log(n), 2) + pol(log(or), 2)),
    h(log(v) ~ pol(log(n), 2) * pol(log(or), 2)) )

  print(ggplot(z, aes(x=error, y=model, shape=type)) + geom_point() +
        facet_wrap(~ stat, scales='free_x') + labs(caption=paste('Model', wg)))
  }
```

The full models (two quadratics) without interaction estimated using robust regression best predicts the log of the standard error of the treatment log odds ratio.  These model have $R^{2} > 0.99$ and less than 0.12 chance of making an error larger than 0.1 on the log scale (which is a fold change error of 0.9).  The median absolute log error is < 0.03 (fold change error of 0.97 in estimating the variance).  Write a function for using these approximations later.


```{r robv}
for(wg in 1 : 2) {
  dat <- subset(w, model == wg)
  approxvar <- MASS::rlm(log(v) ~ pol(log(n), 2) + pol(log(or), 2), data=dat)

  # Derive a function to compute standard error from n and or
  b <- coef(approxvar)
  if(wg == 1) beta.approxvar1 <- b else beta.approxvar2 <- b
  approxse <- function(n, or, n1=NULL, n2=NULL) {
    b <- beta.approxvar
    if(length(n1)) n <- 4. / ((1. / n1) + (1. / n2))
    x <- log(n)
    y <- log(or)
    xb <- b[1] + b[2] * x + b[3] * x ^ 2 + b[4] * y + b[5] * y ^ 2 
    sqrt(exp(xb))
  }
}
```

## Simulation of Frequentist Trials


For each true OR  we simulate 1000 clinical trials.  ORs will vary from 0.2, 0.3, ..., 1.0.  Markov ordinal analysis is done after `r clooks` patients have completed follow-up.   Frequentist power is computed at $\alpha=0.05$.   

For each trial the treatments are randomly assigned with probability 1/2 each.  The high-level `estSeqMarkovOrd` function calls the `simMarkovOrd` function to simulate each trial.   Because the contrast of interest needs to take into account a treatment $\times$ time interaction, the `groupContrast` argument must be specified below.

**Note**: `dosim` calls `estSeqMarkovOrd` which fits data using the `rms` package `lrm` function when the model is a proportional odds model.  Since there is non-proportional odds, the `VGAM` package's `vgam` or `vglm` function is used to fit the partial proportional odds model instead.   When doing Bayesian analysis, the `rmsb` package `blrm` function will fit a constrained partial PO model.

```{r sim, results='asis'}
ppo       <- ~ pol(time, 2)    # partial PO model with PO exception for time
ors       <- c(0.3, 0.4, 0.5, 0.55, 0.6, 0.65, 0.7, 0.8, 1.0)

# Define time to recovery with treatment of death at any time as censored at 28d
tc <- function(y) {
  if(any(y == 2)) return(c(28, 0))   # death
  if(any(y == 0)) return(c(min(which(y == 0)), 1))
	# which works since time=1,2,3,...
  c(length(y), 0)             # censored; length(y) should be 28
}

dosimit  <- function(model, g, formula, file)
  dosim(y=0:2, times=1:28, g=g, initial=c('1'=1.0),  # one initial state, sampled w/prob 1.0
        absorb=c(0, 2),
        intercepts=ints,
        ors=ors, formula=formula,
        ppo=ppo,  groupContrast=if(model == 2) contrvgam,
        timecriterion=tc, coxzph=FALSE, seed=3,
        nsim=1000, looks=looks,
        maxest=2, maxvest=10,
				file=file)

pr('Model 1')
sim1 <- dosimit(1, g1, y ~ pol(time, 2) + group, file='simflooks1.rds')
examSim(sim1, desc='Partial PO Markov model with constant treatment effect')
```

For many of the simulations in this constant treatment effect model the `vgam` function did not converge properly, creating impossibly large maximum likelihood estimates or variance estimates.   In the latter part of the output above these outliers have been removed.  There you can see the extremely good accuracy of the standard error estimates from the maximum likelihood estimation process.

The transition odds ratios are equal to the hazard ratios in the time-to-recovery analyses and at the last data look the power of the transition model treatment test and the Cox proportional hazards model test are the same.  The two $\chi^2$ statistics are virtually the same.

There is approximately 0.94 power to detect An OR or HR of 0.55 and 0.85 power to detect 0.60, all at the planned maximum sample size of 178 for two treatment arms.  As shown earlier, ORs of 0.55 and 0.6 correspond to 3.2 and 2.8 days, respectively, reduction in mean time unrecovered.

Now look at model 2 in which the treatment effect modeled as a quadratic function in time but the true treatment effect is linear up until day 10, then flat.   The transition model treatment contrast is made at 10d using the quadratic time $\times$ treatment interaction. 

```{r sim2,results='asis'}
pr('Model 2')
# beta.approxvar <- beta.approxvar2
sim2 <- dosimit(2, g2, y ~ pol(time, 2) * group, file='simflooks2.rds')
examSim(sim2, desc='Partial PO Markov model with varying treatment effect')
```

Devoting 3 parameters to the treatment effect as done above is very expensive, with major reduction in power when testing the transition odds ratio.  This is entirely consistent with the increased variance of the treatment contrast in the second model.  The next section checks whether this inefficiency carries over to mean time unrecovered.

<a name="cost"></a>

## Number of Treatment Parameters and Efficiency of Reduction in Days Unrecovered

It is possible that as opposed to a treatment contrast of transition odds at a given day, an integrative summary such as mean number of days unrecovered will not inherit the same degree of imprecision / power loss as seen in the variance of the treatment contrast at a specific day when time $\times$ treatment interaction is included in the model.   We explore that here by simulating one dataset with 150 patients in each group, under both model 1 and model 2.  Bayesian Markov proportional odds models are fitted and the posterior distribution of the reduction in mean number of days unrecovered is estimated from 4000 posterior samples.  The data are simulated under a constant treatment effect with OR=0.75 ($\beta=-0.29$) and are analyzed as such in model 1 or using a quadratic time $\times$ treatment interaction in model 2, and requesting mean times in states Y=1-2 over days 1-28.  A constrained partial proportional odds model is used, where the model is in non-PO with regard to time.  With only 3 Y levels, the constrained partial PO model is identical to an unconstrained partial PO model.  The process is repeated for 4 simulated trials.

```{r cost}
require(rmsb)
stanSet()
# Function to compute unconditional probabilities (state occupancy probabilities) for a single
# treatment and model for each posterior draw, then compute all 4000 posterior samples of mean time unrecovered
mtime <- function(fit, tx) {
  u <- soprobMarkovOrdm(fit, data.frame(group=tx, yprev=1), times=1:28, ylevels=0:2, absorb=c(0, 2))
  # Being in state Y=1,2 is equivalent to not being in state 0 (first column of u)
  apply(u, 1, function(z) sum(1. - z[, 1]))
}
ppos <- function(x, eps=0) mean(x > eps)
n <- 300
mt <- NULL
for(isim in 1 : 4) {
  set.seed(10 + isim)
  cat('\nSimulation', isim, '\n\n')
  s1 <- simMarkovOrd(n=n / 2, y=0:2, 1:28, initial=1, X=c(group=1),
     								 absorb=c(0, 2), intercepts=ints, g=g1, parameter=0)
  s2 <- simMarkovOrd(n=n / 2, y=0:2, 1:28, initial=1, X=c(group=2),
      							 absorb=c(0, 2), intercepts=ints, g=g, parameter=log(0.75))
  s       <- rbind(s1, s2)
  s$group <- as.factor(s$group)
  for(wg in 1 : 2) {
    cat(if(wg == 1) 'Constant treatment effect\n\n' else '\nTime-varying treatment effect\n\n')
    form <-  y ~ pol(time, 2) + group
    if(wg == 2) form <- y ~ pol(time, 2) * group
    fi <- paste0('bcost', wg, isim, '.rds')
    b <- blrm(form, ~ pol(time, 2), cppo=function(y) y, data=s, file=fi)
    print(coef(b), digits=3)
    mt <- rbind(mt,  data.frame(model=wg, simulation=isim, delta=mtime(b, '1') - mtime(b, '2')))
  }
}
setDT(mt, key=c('simulation', 'model'))
cat('\nPosterior probabilities of > 0 reduction in days unrecovered:\n\n')
mt[, ppos(delta), by=.(simulation, model)]
cat('\nPosterior probabilities of > 1d reduction in days unrecovered:\n\n')
mt[, ppos(delta, 1), by=.(simulation, model)]

ggplot(mt, aes(x=delta, color=paste('Model', model))) + geom_density() +
  facet_wrap(~ paste('Simulation', simulation), nrow=2) +
  xlab('Reduction in Days Unrecovered') + labs(color='')
```

As judged by the posterior probabilities of any reduction in expected number of days unrecovered and  the width of the posterior distributions for the treatment difference in expected days unrecovered, there is a small to moderate loss in efficiency on the integrative summary when using two more treatment parameters in the transition model.

Now instead of including treatment $\times$ time interactions, assume a constant treatment effect over time but allow for the treatment effect to not be constant over levels of Y, i.e., allow the treatment effect to not satisfy the proportional odds assumption.  We need to define a new `g` function that will simulate such data.  The increment in the treatment log odds ratio for Y=2 is taken as 0.5, meaning that if the treatment is effective for reducing P(Y=1) it is less effective for reducing mortality.  We use the `g3` function defined earlier.

```{r simnonpo}
mt <- NULL
formals(g3)$extra <- c(sparams$extra, 0.5)
## tau=0.5 above; tau ignored when X=1

for(isim in 1 : 4) {
  set.seed(20 + isim)
  cat('\nSimulation', isim, '\n\n')
  s1 <- simMarkovOrd(n=n / 2, y=0:2, 1:28, initial=1, X=c(group=1),
     								 absorb=c(0, 2), intercepts=ints, g=g3, parameter=0)
  s2 <- simMarkovOrd(n=n / 2, y=0:2, 1:28, initial=1, X=c(group=2),
      							 absorb=c(0, 2), intercepts=ints, g=g3, parameter=log(0.75))
  s       <- rbind(s1, s2)
  s$group <- as.factor(s$group)
  for(wg in 1 : 2) {
    cat(if(wg == 1) 'Constant treatment effect with PO\n\n' else '\nCoonstant treatment effect with non-PO\n\n')
    form <-  y ~ pol(time, 2) + group
    ppo  <- ~ pol(time, 2)
    if(wg == 2) ppo <- ~ pol(time, 2) + group
    fi <- paste0('bcostppo', wg, isim, '.rds')
    b <- blrm(form, ppo, cppo=function(y) y, data=s, file=fi)
    print(coef(b), digits=3)
    mt <- rbind(mt,  data.frame(model=wg, simulation=isim, delta=mtime(b, '1') - mtime(b, '2')))
  }
}
setDT(mt, key=c('simulation', 'model'))
cat('\nPosterior probabilities of > 0 reduction in days unrecovered:\n\n')
mt[, ppos(delta), by=.(simulation, model)]
cat('\nPosterior probabilities of > 1d reduction in days unrecovered:\n\n')
mt[, ppos(delta, 1), by=.(simulation, model)]

ggplot(mt, aes(x=delta, color=paste('Model', model))) + geom_density() +
  facet_wrap(~ paste('Simulation', simulation), nrow=2) +
  xlab('Reduction in Days Unrecovered') + labs(color='')
```

# Simulation of Bayesian Sequential Trials

When considering Bayesian operating characteristics, as detailed [here](https://hbiostat.org/bbr/md/alpha.html) it is important to note that the Bayesian approach does not condition on unknowns, so the probability of deciding in favor of efficacy were the true efficacy to be zero is not a probability that makes sense in the Bayesian paradigm.  What is important for a Bayesian analysis to control is the probability of making a wrong decision given the data used **at the decision point**.  This is one minus the probability of efficacy.  This probability can be miscalibrated for only three reasons:

* the data model is incorrect (a problem shared by the frequentist paradigm)
* the study design or measurements are biased (also shared by frequentist)
* the prior is disagreeable (most often caused by an investigator assuming a prior that is too optimistic in favor of a positive treatment benefit)

So the unique aspect to examine in the Bayesian system is the prior distribution.   Once the prior is agreed upon, the accuracy of posterior probabilities flows automatically, given the first two bullet points above are not problematic.

For those who insist on applying frequentist concepts to interpretation of Bayesian results, Bayesian power and the Bayesian equivalent of type I assertion probability $\alpha$ under the null (OR=1.0) are given below.

Since we are using the normal approximation to the data likelihood, the time-consuming part of the simulations has already be done.  What remains is to translate the log odds ratio estimates their estimated variances to posterior distributions.  We use flat priors for all but the treatment effect, and apply a skeptical normal prior distribution for the treatment effect.  Then the posterior distribution for the treatment effect is also normal.  The skeptical prior is set up as follows.  We assume it is equally likely that the treatment causes harm as that it benefits the patient.  Thus the mean of the normal prior distribution for the treatment effect log(OR) is zero.  The standard deviation is chosen so that P(OR < 0.25) = P(OR > 4) = 0.05.  Three other priors are used for assessing efficacy: a more skeptical one, an even more skeptical one, and one that is nearly flat.

The data model assumes a constant treatment effect on the transition OR scale.

Posterior probabilities will be computed for the following assertions:

* Efficacy: P(OR < 1)
* Efficacy: P(OR < 1) with a target of 0.975 instead of 0.95 (all other efficacy assessments target 0.95)
* Efficacy, more skeptical prior: P(OR < 0.25) = P(OR > 4) = 0.01
* Efficacy, extremely skeptical prior: P(OR < 0.25) = P(OR > 4) = 0.001
* Efficacy, nearly flat prior: P(OR < 1) with prior P(OR < 0.25) = P(OR > 4) = 0.4.
* Inefficacy/harm: P(OR > 1)


Bayesian power is taken to be the probability that the posterior probability of a beneficial effect exceeds 0.95 at the end of the study.  We also compute the probability of finding a greater than 0.9 probability of a harmful effect.

```{r bayespost}
asserts <-
  list(list('Efficacy',                       '<', 0, cutprior=log(4), tailprob=0.05),
       list('Efficacy, target 0.975',         '<', 0, cutprior=log(4), tailprob=0.01),
       list('Efficacy, more skeptical prior', '<', 0, cutprior=log(4), tailprob=0.01),
       list('Efficacy, extremely skeptical prior', '<', 0, cutprior=log(4), tailprob=0.001),
       list('Efficacy, flat prior',           '<', 0, cutprior=log(4), tailprob=0.4),
       list('Inefficacy/harm',                '>', 0, cutprior=log(4), tailprob=0.05))
  
s <- gbayesSeqSim(sim1$sim, asserts=asserts)
setDT(s, key=c('OR', 'sim', 'look'))
s[, c('parameter', 'loghr', 'lrchisq') := NULL]
head(s)
attr(s, 'asserts')
alabels <- attr(s, 'alabels')  # named vector to map p1 p2 to labels
```

Compute the Bayesian power---the probability of hitting assertion-specific targets at the planned study end.

```{r hits,fig.cap='Probability of hitting posterior probability targets for various assertions, based on the posterior distribution of the odds ratio.  Posterior probability targets could also be constructed for state occupancy probabilities.'}
# Reshape results into taller and thinner data table so can plot over multiple assertions

ps <- names(alabels)
m  <- melt(s,
           measure.vars=list(ps, paste0('mean', 1:5), paste0('sd', 1:5)),
           variable.name='assert', value.name=c('p', 'mean', 'sd'))
m[, assert := alabels[assert]]
setkey(m, OR, assert, sim, look)
head(m)

# Re-check frequentist power
z <- m[assert == 'Efficacy' & look == max(looks), ]
pr('Frequentist power',
   z[, .(Power=mean(abs(est / sqrt(vest)) > 1.959964, na.rm=TRUE)), by=OR])
# Check Bayesian power if only looked at the end
pr('Bayesian power with one data look at n=178',
   z[, .(Power=mean(p > 0.95, na.rm=TRUE)), by=OR])

# Define targets
ptarget <- c('Efficacy'                            = 0.95,
             'Efficacy, target 0.975'              = 0.975,
             'Efficacy, more skeptical prior'      = 0.95,
             'Efficacy, extremely skeptical prior' = 0.95,
             'Efficacy, flat prior'                = 0.95,
             'Inefficacy/harm'                     = 0.9 )

m[, ptarget := ptarget[assert]]	  # spreads targets to all rows

# Show the posterior probability of efficacy paths for 12 of the simulated RCTs
w <- m[assert == 'Efficacy' & sim < '1 20', ]
setkey(w, OR, sim, look)
ggplot(w, aes(look, p, color=factor(sim))) + geom_line() +
  facet_wrap(~ OR) + xlab('Sample Size at Data Look') +
  ylab('Posterior Probability of Efficacy') +
  labs(color='', title='Representative posterior probability of efficacy paths',
       subtitle='12 simulated RCTs, one panel per true OR') +
  guides(color=FALSE)

# Target is hit if the posterior probability of the assertion ever reaches
# the probability target.   This happens if and only if the maximum
# post. prob. exceeds the target, so compute pmax for each combination of
# OR, assert, and sim

pm <- m[, .(pmax = max(p, na.rm=TRUE), ptarget=ptarget[1]), by=.(OR, assert, sim)]
# Compute mean maximum post prob
pmm <- pm[, .(meanpp = mean(pmax, na.rm=TRUE)), by=.(OR, assert)]
ggplot(pmm, aes(x=OR, y=meanpp, color=assert)) + geom_point() +
  labs(color='', title='Mean Maximum Posterior Probability over Data Looks',
       subtitle='Evidence for harm is never great since no OR > 1 were tested')

# hit = 0/1 indicator if hitting target at the single fixed sample size
u <- pm[, .(Power = mean(pmax > ptarget, na.rm=TRUE)), by=.(OR, assert)]
z <- u
z$txt <- with(z, paste0('OR:', OR, '<br>', assert, '<br>Hit probability:', round(Power, 3)))
ggp(ggplot(z, aes(x=OR, y=Power, color=assert, label=txt)) +
  geom_line() +
  xlab('OR') + ylab('Proportion Hitting Posterior Probability Target') +
  guides(color=guide_legend(title='Assertion'))) 
```

From the table below, Bayesian power based on multiple data looks is 0.960 if the true OR=0.6 and the original prior distribution is used.   Under the more skeptical prior power drops to 0.955.  The Bayesian equivalent of the type I assertion probability $\alpha$ is found by viewing the row with OR=1.0, where the proportion of simulated clinical trials that would hit the 0.95 posterior probability efficacy target at any look is 0.143, and is 0.123 for the more skeptical prior.  The sacrifice in Bayesian power from moving to the more skeptical prior distribution is a difference of 0.005 at OR=0.6.  The loss is a little greater at OR=0.8.  Instead of using a posterior probability cutoff of 0.95, when 0.975 is used in conjunction with the original skeptical prior, the probability if ever hitting the efficacy target when OR=1.0 is 0.062.

```{r powertable}
uu <- dcast(u, OR ~ assert, value.var='Power')
knitr::kable(uu, digits=3, caption='Bayesian power: probability of hitting the target at any look')
```

# Futility Analysis

Futility is based on the probability that, given the current posterior probability $P_c$ at one of the first data looks, the posterior probability of efficacy will never be $\geq 0.95$, written symbolically as $P(\max(P(OR < 1) | P_{c})) < 0.95)$ where the maximum is over the remaining data looks.  This failure probability will be estimated as a function of $P_c$ by simulation, and is not conditional on a value of the OR.

More details about futility calculations may be found [here](https://hbiostat.org/R/Hmisc/gbayesSeqSim.html).

For looks at n=`r clooks1` save the posterior probabilities of efficacy and use these to predict, using a binary logistic model, the probability that the posterior probability is $\geq$ 0.95.

```{r fut}
w <- s[, .(OR, sim, look, p1)]
# Save final post prob and put on each record
w[, p := p1[.N], by=.(OR, sim)]
w <- w[look < max(looks),]
w[, n := factor(look)]
dd <- datadist(w); options(datadist='dd')
cqnorm <- function(x) qnorm(pmax(pmin(x, 1. - 1e-5), 1e-5))
f <- lrm(p >= 0.95 ~ rcs(cqnorm(p1), 4) * n, data=w, maxit=30)
ggp(ggplot(Predict(f, p1, n, fun=plogis, conf.int=FALSE), colorscale=function(...) {}) +
  scale_x_continuous(breaks=seq(0, 1, by=0.1), minor_breaks=seq(0, 1, by=0.05)) +
  scale_y_continuous(breaks=seq(0, 1, by=0.1), minor_breaks=seq(0, 1, by=0.05)) +
  xlab('Early Posterior Probability') +
  ylab('P(Final posterior probability >= 0.95)'))
```

For there to be better than a 50/50 chance of ultimate success the posterior probability must exceed 0.4 at the first look (n=20 for two treatment arms) or 0.86 at n=100.

See what proportion of trials, separately for each true unknown OR, would have stopped early for futility, here meaning the probability that the ultimate posterior probability exceeds 0.95 being less than 0.05.

```{r fut2}
# Create an R function that computes the prob. that the last post. prob. exceeds 0.95
success <- Function(f)
failure  <- function(pp, n) 1. - plogis(success(pp, as.character(n)))
failure(0.4, 20); failure(0.86, 100)
# Put failure probability on each look
w[, fail := failure(p1, look)]
z <- w[, .(failed=any(fail > 0.95)), by=.(OR, sim)]
pr('Proportion of trials stopped early for futility',
   z[, .(propfailed = mean(failed, na.rm=TRUE)), by=OR])
```

There is a substantial (0.48) risk of stopping early for futility were the true OR to be $\geq 0.8$.  Under the null, thankfully 0.90 of trials would stop early for futility, here futility meaning low chance of ever declaring efficacy.

Were the Bayesian operating characteristics to be re-run with stopping early when the ultimate success probability is $< 0.05$, the $\alpha$ probability would be lower, but the Bayesian power would also be noticeably lower were OR to be $\geq 0.65$ (where 0.13 of trials would have stopped early without claiming efficacy).

# More Information

* Full R markdown script, `Hmisc` and service functions, and `pdf` version of this report: [hbiostat.org/R/Hmisc/markov](https://hbiostat.org/R/Hmisc/markov)
* COVID-19 statistical resources: [hbiostat.org/proj/covid19](https://hbiostat.org/proj/covid19) including detailed analyses of the VIOLET 2 and ORCHID studies plus an examination of the handling of irregular measurement times in a Markov model
* Bayesian design and analysis resources: [hbiostat.org/bayes](https://hbiostat.org/bayes)
* Markov modeling references: [hbiostat.org/bib/markov.html](https://hbiostat.org/bib/markov.html)
* Longitudinal ordinal analysis references: [hbiostat.org/bib/ordSerial.html](https://hbiostat.org/bib/ordSerial.html)
* General references on longitudinal data analysis: [hbiostat.org/bib/serial.html](https://hbiostat.org/bib/ordSerial.html)
* Introduction to the proportional odds model: [hbiostat.org/bbr/md](https://hbiostat.org/bbr/md) chapter 7
* More about the proportional odds model: [hbiostat.org/rms](https://hbiostat.org/rms)

# Computing Environment

`r markupSpecs$html$session()`

```{r quit,echo=FALSE}
# knitr::knit_exit()
```
