1 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. 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.

As shown below, 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.

2 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. As detailed here, 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 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) partial proportional odds model.

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 30, 54, 78, 102, 126, 150, 174, 198, 222, 246, and 267 (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 20, 36, 52, 68, 84, 100, 116, 132, 148, 164, and 178 (final look at planned maximum sample size).

3 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
}

4 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, 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.

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.

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)

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

[1] 0

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

[1] 0

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)

[1] 0.704

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.

# 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.

5 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.

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.

6 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.

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 0.702 to 0.855.

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

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.

8 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.

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

Constant treatment effect

beta = -0.638 

Time-varying treatment effect

$Contrast
   group2 
-1.129932 

$SE
[1] 0.00553126

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

Median model 2 : model 1 ratio of variances of treatment effect NA

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

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

Model 1 

Model 2 

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

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

Model 1 

sim1 <- dosimit(1, g1, y ~ pol(time, 2) + group, file='simflooks1.rds')

Failed simulations 

                                                                                               failure
1:                                                                                        |contrast|>2
2:                                        Error in lm.fit(X.vlm, y = z.vlm, ...) : NA/NaN/Inf in 'y'\n
3:              Error in tapplymat1(ccump, "diff") : \n  NA/NaN/Inf in foreign function call (arg 1)\n
4: Simulated data for simulation with sample size 20 has 2 distinct y values instead of the required 3
5: Simulated data for simulation with sample size 36 has 2 distinct y values instead of the required 3
6: Simulated data for simulation with sample size 52 has 2 distinct y values instead of the required 3
    Freq
1:   131
2:  3485
3: 16997
4:   168
5:    10
6:     1

Median of 1000 logistic regression coefficients for each OR 

     (Intercept):1 (Intercept):2 pol(time, 2)time:1 pol(time, 2)time:2
0.3         5.0555       -3.7088            -0.1263            -0.2686
0.4         4.7127       -3.6245            -0.0437            -0.3069
0.5         4.7273       -3.8868            -0.0757            -0.2698
0.55        4.8425       -3.6212            -0.0772            -0.2361
0.6         5.0698       -3.8890            -0.1264            -0.2352
0.65        4.8132       -3.8957            -0.0644            -0.1989
0.7         4.7954       -3.8270            -0.0586            -0.2001
0.8         4.7900       -3.8663            -0.0979            -0.2190
1           4.9208       -4.0593            -0.1362            -0.2037
     pol(time, 2)time^2:1 pol(time, 2)time^2:2  group2
0.3               -0.0047               0.0176 -0.8863
0.4               -0.0088               0.0180 -0.7568
0.5               -0.0064               0.0162 -0.5733
0.55              -0.0060               0.0142 -0.5798
0.6               -0.0034               0.0160 -0.4835
0.65              -0.0059               0.0130 -0.4316
0.7               -0.0070               0.0124 -0.2652
0.8               -0.0042               0.0130 -0.2146
1                 -0.0026               0.0132  0.0322

Standard deviation of simulated group effects, scaled MAD, and square root of median of estimated variances 

      OR look    SD   MAD SDest
 1: 0.30   20 0.397 0.188 0.513
 2: 0.30   36 0.282 0.130 0.244
 3: 0.30   52 0.239 0.120 0.186
 4: 0.30   68 0.208 0.104 0.156
 5: 0.30   84 0.210 0.089 0.138
 6: 0.30  100 0.204 0.080 0.123
 7: 0.30  116 0.206 0.086 0.006
 8: 0.30  132 0.196 0.085 0.001
 9: 0.30  148 0.192 0.087 0.000
10: 0.30  164 0.185 0.087 0.000
11: 0.30  178 0.180 0.089 0.000
12: 0.40   20 0.345 0.209 0.339
13: 0.40   36 0.163 0.176 0.004
14: 0.40   52 0.127 0.174 0.000
15: 0.40   68 0.123 0.174 0.000
16: 0.40   84 0.118 0.174 0.000
17: 0.40  100 0.117 0.174 0.000
18: 0.40  116 0.119 0.174 0.000
19: 0.40  132 0.117 0.174 0.000
20: 0.40  148 0.117 0.174 0.000
21: 0.40  164 0.116 0.174 0.000
22: 0.40  178 0.119 0.174 0.000
23: 0.50   20 0.413 0.215 0.389
24: 0.50   36 0.235 0.161 0.004
25: 0.50   52 0.182 0.181 0.000
26: 0.50   68 0.171 0.206 0.000
27: 0.50   84 0.175 0.205 0.000
28: 0.50  100 0.174 0.205 0.000
29: 0.50  116 0.172 0.203 0.000
30: 0.50  132 0.173 0.206 0.000
31: 0.50  148 0.171 0.206 0.000
32: 0.50  164 0.172 0.206 0.000
33: 0.50  178 0.171 0.206 0.000
34: 0.55   20 0.463 0.362 0.518
35: 0.55   36 0.292 0.168 0.010
36: 0.55   52 0.214 0.164 0.001
37: 0.55   68 0.184 0.151 0.000
38: 0.55   84 0.171 0.142 0.000
39: 0.55  100 0.159 0.143 0.000
40: 0.55  116 0.155 0.141 0.000
41: 0.55  132 0.148 0.140 0.000
42: 0.55  148 0.145 0.140 0.000
43: 0.55  164 0.143 0.140 0.000
44: 0.55  178 0.141 0.140 0.000
45: 0.60   20 0.494 0.510 0.519
46: 0.60   36 0.314 0.240 0.368
47: 0.60   52 0.265 0.232 0.005
48: 0.60   68 0.237 0.230 0.004
49: 0.60   84 0.220 0.205 0.003
50: 0.60  100 0.205 0.204 0.002
51: 0.60  116 0.195 0.182 0.000
52: 0.60  132 0.188 0.168 0.000
53: 0.60  148 0.184 0.156 0.000
54: 0.60  164 0.179 0.158 0.000
55: 0.60  178 0.177 0.145 0.000
56: 0.65   20 0.394 0.205 0.319
57: 0.65   36 0.213 0.095 0.000
58: 0.65   52 0.165 0.096 0.000
59: 0.65   68 0.141 0.095 0.000
60: 0.65   84 0.130 0.095 0.000
61: 0.65  100 0.121 0.076 0.000
62: 0.65  116 0.116 0.052 0.000
63: 0.65  132 0.114 0.054 0.000
64: 0.65  148 0.112 0.052 0.000
65: 0.65  164 0.110 0.051 0.000
66: 0.65  178 0.110 0.051 0.000
67: 0.70   20 0.404 0.256 0.017
68: 0.70   36 0.190 0.165 0.000
69: 0.70   52 0.144 0.163 0.000
70: 0.70   68 0.134 0.163 0.000
71: 0.70   84 0.133 0.163 0.000
72: 0.70  100 0.135 0.163 0.000
73: 0.70  116 0.132 0.163 0.000
74: 0.70  132 0.131 0.163 0.000
75: 0.70  148 0.131 0.163 0.000
76: 0.70  164 0.131 0.163 0.000
77: 0.70  178 0.130 0.163 0.000
78: 0.80   20 0.511 0.323 0.507
79: 0.80   36 0.296 0.211 0.011
80: 0.80   52 0.219 0.181 0.006
81: 0.80   68 0.189 0.181 0.005
82: 0.80   84 0.179 0.181 0.005
83: 0.80  100 0.173 0.181 0.004
84: 0.80  116 0.163 0.181 0.004
85: 0.80  132 0.162 0.181 0.004
86: 0.80  148 0.158 0.178 0.003
87: 0.80  164 0.155 0.181 0.003
88: 0.80  178 0.151 0.165 0.003
89: 1.00   20 0.523 0.441 0.516
90: 1.00   36 0.330 0.245 0.365
91: 1.00   52 0.252 0.227 0.006
92: 1.00   68 0.226 0.226 0.002
93: 1.00   84 0.203 0.224 0.001
94: 1.00  100 0.187 0.203 0.001
95: 1.00  116 0.177 0.202 0.001
96: 1.00  132 0.166 0.202 0.000
97: 1.00  148 0.163 0.202 0.000
98: 1.00  164 0.157 0.202 0.000
99: 1.00  178 0.155 0.199 0.000
      OR look    SD   MAD SDest

examSim(sim1, desc='Partial PO Markov model with constant treatment effect')

Number of successful simulations out of 1000 

      OR look   n
 1: 0.30   20 551
 2: 0.30   36 646
 3: 0.30   52 702
 4: 0.30   68 741
 5: 0.30   84 754
 6: 0.30  100 774
 7: 0.30  116 791
 8: 0.30  132 800
 9: 0.30  148 815
10: 0.30  164 822
11: 0.30  178 831
12: 0.40   20 471
13: 0.40   36 545
14: 0.40   52 636
15: 0.40   68 699
16: 0.40   84 761
17: 0.40  100 817
18: 0.40  116 848
19: 0.40  132 872
20: 0.40  148 886
21: 0.40  164 902
22: 0.40  178 915
23: 0.50   20 399
24: 0.50   36 472
25: 0.50   52 525
26: 0.50   68 585
27: 0.50   84 646
28: 0.50  100 688
29: 0.50  116 734
30: 0.50  132 759
31: 0.50  148 779
32: 0.50  164 796
33: 0.50  178 808
34: 0.55   20 528
35: 0.55   36 635
36: 0.55   52 736
37: 0.55   68 802
38: 0.55   84 856
39: 0.55  100 895
40: 0.55  116 915
41: 0.55  132 931
42: 0.55  148 945
43: 0.55  164 957
44: 0.55  178 962
45: 0.60   20 569
46: 0.60   36 633
47: 0.60   52 698
48: 0.60   68 770
49: 0.60   84 798
50: 0.60  100 835
51: 0.60  116 860
52: 0.60  132 871
53: 0.60  148 881
54: 0.60  164 884
55: 0.60  178 888
56: 0.65   20 460
57: 0.65   36 586
58: 0.65   52 688
59: 0.65   68 760
60: 0.65   84 816
61: 0.65  100 840
62: 0.65  116 863
63: 0.65  132 876
64: 0.65  148 887
65: 0.65  164 891
66: 0.65  178 893
67: 0.70   20 501
68: 0.70   36 614
69: 0.70   52 742
70: 0.70   68 822
71: 0.70   84 879
72: 0.70  100 915
73: 0.70  116 943
74: 0.70  132 961
75: 0.70  148 974
76: 0.70  164 981
77: 0.70  178 984
78: 0.80   20 542
79: 0.80   36 605
80: 0.80   52 708
81: 0.80   68 777
82: 0.80   84 818
83: 0.80  100 847
84: 0.80  116 869
85: 0.80  132 888
86: 0.80  148 894
87: 0.80  164 901
88: 0.80  178 905
89: 1.00   20 629
90: 1.00   36 745
91: 1.00   52 841
92: 1.00   68 910
93: 1.00   84 943
94: 1.00  100 961
95: 1.00  116 975
96: 1.00  132 984
97: 1.00  148 986
98: 1.00  164 992
99: 1.00  178 993
      OR look   n

Comparison of OR and HR 

     OR   HR
1: 0.30 0.34
2: 0.40 0.43
3: 0.50 0.52
4: 0.55 0.56
5: 0.60 0.60
6: 0.65 0.65
7: 0.70 0.70
8: 0.80 0.80
9: 1.00 1.00

Power of Cox test at last look for time until Y=1 

     OR power
1: 0.30  1.00
2: 0.40  1.00
3: 0.50  0.96
4: 0.55  0.92
5: 0.60  0.84
6: 0.65  0.69
7: 0.70  0.53
8: 0.80  0.26
9: 1.00  0.04

Power of Markov proportional odds model test 

      OR look power
 1: 0.30   20  0.70
 2: 0.30   36  0.93
 3: 0.30   52  0.99
 4: 0.30   68  1.00
 5: 0.30   84  1.00
 6: 0.30  100  1.00
 7: 0.30  116  1.00
 8: 0.30  132  1.00
 9: 0.30  148  1.00
10: 0.30  164  1.00
11: 0.30  178  1.00
12: 0.40   20  0.70
13: 0.40   36  0.93
14: 0.40   52  0.98
15: 0.40   68  0.99
16: 0.40   84  1.00
17: 0.40  100  1.00
18: 0.40  116  1.00
19: 0.40  132  1.00
20: 0.40  148  1.00
21: 0.40  164  1.00
22: 0.40  178  1.00
23: 0.50   20  0.53
24: 0.50   36  0.88
25: 0.50   52  0.97
26: 0.50   68  0.99
27: 0.50   84  0.99
28: 0.50  100  1.00
29: 0.50  116  1.00
30: 0.50  132  1.00
31: 0.50  148  1.00
32: 0.50  164  1.00
33: 0.50  178  1.00
34: 0.55   20  0.44
35: 0.55   36  0.71
36: 0.55   52  0.84
37: 0.55   68  0.90
38: 0.55   84  0.93
39: 0.55  100  0.95
40: 0.55  116  0.96
41: 0.55  132  0.97
42: 0.55  148  0.98
43: 0.55  164  0.98
44: 0.55  178  0.99
45: 0.60   20  0.36
46: 0.60   36  0.60
47: 0.60   52  0.72
48: 0.60   68  0.79
49: 0.60   84  0.84
50: 0.60  100  0.87
51: 0.60  116  0.90
52: 0.60  132  0.92
53: 0.60  148  0.94
54: 0.60  164  0.95
55: 0.60  178  0.95
56: 0.65   20  0.54
57: 0.65   36  0.81
58: 0.65   52  0.89
59: 0.65   68  0.92
60: 0.65   84  0.93
61: 0.65  100  0.94
62: 0.65  116  0.95
63: 0.65  132  0.96
64: 0.65  148  0.96
65: 0.65  164  0.97
66: 0.65  178  0.97
67: 0.70   20  0.56
68: 0.70   36  0.87
69: 0.70   52  0.96
70: 0.70   68  0.98
71: 0.70   84  0.98
72: 0.70  100  0.99
73: 0.70  116  0.99
74: 0.70  132  1.00
75: 0.70  148  1.00
76: 0.70  164  1.00
77: 0.70  178  0.99
78: 0.80   20  0.43
79: 0.80   36  0.64
80: 0.80   52  0.72
81: 0.80   68  0.76
82: 0.80   84  0.78
83: 0.80  100  0.80
84: 0.80  116  0.81
85: 0.80  132  0.82
86: 0.80  148  0.83
87: 0.80  164  0.84
88: 0.80  178  0.84
89: 1.00   20  0.35
90: 1.00   36  0.52
91: 1.00   52  0.59
92: 1.00   68  0.63
93: 1.00   84  0.64
94: 1.00  100  0.64
95: 1.00  116  0.65
96: 1.00  132  0.65
97: 1.00  148  0.65
98: 1.00  164  0.65
99: 1.00  178  0.65
      OR look power

Spearman rho correlation between Markov and Cox model chi-square: 0.24