---
title: Simulating Operating Characteristics of Longitudinal Markov Ordinal Outcome
  Trials
author: |
  | Frank Harrell
  | VICTR Research Methods Program
  | Department of Biostatistics
  | Vanderbilt University School of Medicine
  | CONNECTS VUMC Study Design and Analysis Core 5
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
  html_document:
    toc_depth: '3'
    df_print: paged
  pdf_document:
    latex_engine: lualatex
header-includes: \usepackage{longtable}
link-citations: yes
description: Markov Simulation
---
<!-- To compile: Rstudio knit  -->
<style>
.table {
	width: 50%;
}
p.caption {
  font-size: 0.8em;
  color: DarkBlue;
}
</style>

```{r setup,echo=FALSE,include=FALSE}
outfmt <- c('html', 'pdf')[1]
markup <- switch(outfmt, html='html', pdf='latex')
require(rms)          # automatically engages rms(Hmisc)
source('simMarkovOrd.r')
require(data.table)   # uses to compute various summary measures
fdev <- switch(outfmt, html='png', pdf='pdf')
knitrSet(lang='markdown', w=7, h=7, dev=fdev, fig.path=paste0(fdev, '/sim-'))
options(prType=markup)

simdone        <- TRUE
simdone.single <- TRUE

pr <- function(x, obj=NULL, inline=NULL) {
  cat(x, inline, '\n\n', sep='')
	if(length(obj)) print(obj)
	invisible()
	}
```

# Introduction

As opposed to time-to-event analysis or comparing outcomes at a single follow-up time, longitudinal ordinal analysis makes optimum use of available data from multiple time periods, resulting in greater Bayesian and frequentist statistical power or smaller sample sizes to achieve the same power as information-losing methods.  This report provides a template for simulating operating characteristics for treatment comparisons in a longitudinal study with an ordinal outcome Y where the serial correlation is modeled using a first-order Markov process.  This within-patient correlation model is specified by conditioning on the ordinal outcome at the previous time point just as though it were any covariate, with one exception.  Since the time gaps between measurements (which are made on days 1, 3, 7, 14, 28 in these examples) are not constant, an interaction term is added to the model so that the influence of the previous state wanes as the time gap increases.

Results for both a frequentist analysis and a Bayesian analysis are provided.  Only one data look is taken, at the final sample size.  For the Bayesian calculations, the simulations presented here assume that a normal distribution is an adequate approximation to the distribution of the treatment effect estimate, here a treatment B : treatment A transition log odds ratio.  The data model is a proportional odds ordinal logistic model analyzed with the R `rms` package `lrm` function.  Simulations are done by the `Hmisc` package `simMarkovOrd` and `estSeqMarkovOrd` [functions](https://github.com/harrelfe/Hmisc/blob/master/R/simMarkovOrd.r).

The Markov proportional odds ordinal logistic model used here is stated in terms of the probability of being in state $y$ or worse given baseline covariates (including treatment) and the previous state.  This state transition model has the following form: $\Pr(Y(t) \geq y | X, Y(t-\delta)) = \mathrm{expit}(\alpha_{y} + X\beta + f(t, \delta, Y(t-\delta)))$ where $t-\delta$ is the most recent measurement time before current time $t$, $f$ is a linear model involving a number of regression coefficients that are fitted in addition to $\alpha$ and $\beta$, and expit is the logistic cumulative distribution function (inverse logit).  When $t$ is the first post-time-zero measurement, $Y(t-\delta)$ is the baseline state.

State transition probabilities are easy to understand, but for study planning and final interpretation we are more used to seeing state occupancy probabilities (SOP).  When a state is an absorbing state (e.g., death), the SOP is the cumulative incidence of that state.  In general, a SOP is the probability of being in a state $Y=y$ **at** a given time $t$, whether or not the patient was in a certain state previously.  To compute SOPs we uncondition on previous states so that the result is conditional only on the pre-study state and baseline covariates $X$.  In this context, the `Hmisc` function `soprobMarkovOrd` is exemplified.  This function computes exact SOPs from the simulation model.  The `Hmisc` `intMarkovOrd` function is used to compute the proportional odds model intercepts that satisfy specified SOPs for Y=1, 2, 3, 4 at day 28.  These intercepts are used in the simulations.

The new functions are in `Hmisc` version 4.4-3.  Until this is available from CRAN, source, Linux binary, and Windows binary versions are available at [hbiostat.org/R/Hmisc](https://hbiostat.org/R/Hmisc).

As implemented in the R `Hmisc` package [`gbayes`](https://www.rdocumentation.org/packages/Hmisc/versions/4.4-1/topics/gbayes) function, the Bayesian posterior distribution used to approximate (thus avoiding another simulation loop for MCMC posterior draws) the real posterior distribution is normal, making for very quick approximate posterior probability calculations when the prior distribution used for the log OR is also Gaussian as used here.  Three priors are used:

* a skeptical prior for assessing evidence for efficacy
* a flat prior for assessing evidence for efficacy
* a flat prior for studying posterior probabilities of inefficacy/harm
* an optimistic prior for inefficacy/harm

See [hbiostat.org/proj/covid19/bayesplan.html](https://hbiostat.org/proj/covid19/bayesplan.html) and [hbiostat.org/stat/irreg.html](https://hbiostat.org/stat/irreg.html) for more simulations and graphical presentations of them.

The `estSeqMarkovOrd` function also optionally provides, at the last data look only, the Cox proportional hazards $\chi^2$ statistic for treatment for each simulated clinical trial, and the $\chi^2$ statistic for testing proportional hazards in this unadjusted Cox model.  The Cox test for treatment is done after the simulated serial ordinal responses for a patient are summarized with the time until achieving Y=1, with this time right censored at 28 if the patient never achieved Y=1 within 28d.  The purpose of this additional simulation is to compare the frequentist power of a two-sample comparison of time until Y=1 with the frequentist power of the transition odds ratio at day 28 in the Markov proportional odds model, and to gauge the extent to which the time-to-event variable induced by the Markov model operates in proportional hazards across the two treatment groups.

Finally, this report presents a simulation of frequentist power for comparing ordinal outcomes at a single day of follow-up, for a single odds ratio.

# Simulation Parameters

The ordinal outcome variable is taken to have values Y=1,2,3,4 measured on days 1, 3, 7, 14, 28.

The treatment effect on transition probabilities, on the log odds ratio scale, is taken to be linear in time, i.e., the treatment effect is zero on day 1 and linearly increases to its maximum on day 28.  This maximum value is used in the graphical output and in specifying the underlying true treatment effect.
Since a large number of true ORs are used in the simulation, we  simulate 1000 clinical trials per OR.  ORs will vary from benefit to harm over this sequence: 0.4, 0.5, 0.6, ..., 1.0, 1.25.  The sample size is 600 patients.  Frequentist power is computed at $\alpha=0.05$.   Posterior probabilities will be computed for the following assertion and prior distribution combinations:

* Efficacy: P(OR < 1) with a skeptical prior P(OR < 0.5) = 0.05
* Efficacy: P(OR < 1) with a flat prior (log(OR) has mean 0 and SD 100)
* Inefficacy/harm: P(OR > 1) with a flat prior (log(OR) has mean 0 and SD 100)
* Inefficacy/harm: P(OR > 1) with an optimistic prior (log(OR) has mean log(0.85) = -0.1625 and SD of 0.5)

# Simulation Model Specification

## Simulation Model Components

This report addresses a 4-level ordinal outcome Y with levels $y=1,2,3,4$.  These levels are numeric for simplicity but the numeric values are never used.  We might just as well have coded $y=a,b,c,d$.  Throughout we use measurement times $t=1,2,3,7,14,28$.

Careful specification of longitudinal simulation models results in realistic simulated data that are useful for computing sample sizes and operating characteristics.  For Markov transition models there are several components to consider:

* Aspects which are in common with all longitudinal models: choice of measurement times, treatment allocation ratio, and baseline covariate distributions
* The distribution of a special baseline covariate: the initial patient state $Y(0)$
* The intercepts $\alpha_{y}$ in the proportional odds model corresponding to $y=2,3,4$ (there is no intercept for $\Pr(Y \geq 1)$ since this probability is 1.
* The treatment effect $\beta$ which is under our control, and is the log odds ratio for state transition probabilities.  Our simulation model assumes that the treatment effect (log odds ratio for group 2 : group 1) is zero at the first follow-up measurement $t=1$ and that it linearly increases to the last day of follow-up, day 28.  So the treatment effect is $\beta \frac{t - 1}{27}$ for simulating data.  For model fitting with real data we would estimate two parameters: the treatment main effect and the treatment $\times t$ interaction effect.  Were the time $\times$ treatment interaction not linear there would be more terms.  Simulations below explore the power loss from adding an unnecessary quadratic term to this part of the model.
* Imposition of some states being absorbing so that the patient remains in that state forever with probability 1.  Our examples take $y=4$ to be an absorbing state (e.g., death).
* The effects of the previous state on the log odds that the current state is $\geq y$; this is how within-subject dependence is modeled.  There are three possible previous states: $y=1,2,3$ since once $y=4$ is reached there are no more records for the patient.  On the log odds scale we write this part of the model as $\tau_{1} [y' = 1] + \tau_{2} [y' = 2] + \tau_{3} [y' = 3]$, where the $\tau$s are to be determined and $[c]$ denotes 1 if $c$ holds, 0 otherwise.
* When the spacing being measurement times is not constant, we need to model how the influence of the previous state may wane as the time gap from that measurement increases.  We choose to multiply this discounting function by the previous state influence function in the last item, and take the discounting function to be $\exp(-\gamma (\delta - 1))$ where the decay constant $\gamma$ is to be determined and $\delta$ is the observed time gap from the previous obervation.  This formulation assumes that the minimum time gap is one unit, and when $\delta=1$ the previous state gets full weight.
* A general time trend indicating that the majority of patients tend to get more well regardless of treatment or initial state.  This trend is assumed to be lienar in our simulations so is specified as $\kappa (t - 1)$.

The only covariate $X$ in our model is treatment.  Putting all this together, our Markov transition model is

$$\Pr(Y(t) \geq y | Y(t - \delta), X) = \mathrm{expit}(\alpha_{y} + \beta [X=2] + \kappa (t -1) + \exp(-gamma (\delta - 1)) (\tau_{1}[Y(t-\delta) = 1] + \tau_{2}[Y(t-\delta) = 2] + \tau_3{Y(t-\delta) = 3])$$
where $\mathrm{expit}$ is the inverse logit transformation or $\frac{1}{1 + \exp(-x}$.

## Setting Simulation Model Parameter Values

The next step is to specify characteristics to achieve with the simulation model.  We start the process by setting SOPs at $t=1$ to constants when a patient is in treatment group 1 and starts at a given initial state.  We use as our reference an initial state of $y=2$ which is taken to be the most common state at baseline.  For group 1 we ignore $\beta$.  The model reduces to $\Pr(Y(1) \geq y | Y(0), X) = \mathrm{expit}(\alpha_{y} + \tau_{2}$.  We have 4 parameters to solve for given three constraints corresponding to required SOPs for $y=1,2,3,4$ (three constraints since the sum of the four SOPs must be 1).  If we take $\tau_{2}=0$ we can solve for $\alpha$ easily.  Suppose that we require for group 1 that the SOPs at $t=1$ be $[0.05, 0.70, 0.24, 0.01]$ corresponding to $y=1,2,3,4$.  The cumulative probabilities for $Y\geq y$ for $y=2,3,4$ are 0.95, 0.25, 0.01 and their logits are `r round(qlogis(c(0.95, 0.25, 0.01)), 2)`.  We will take these as are starting values for $\alpha$ when $Y(0)=2$.

We must constrain the model using another time point in order to solve for the state dependencies.  We take that second time to be $t=28$, i.e., the last planned follow-up time.  Aside from the three $\alpha$s, these parameters need to be solved for: $kappa, \gamma, \tau_{1}, \tau_{2}, \tau_{3}$.  If again we force $\tau_{2}=0$ we still have 4 parameters for 3 constraints.  We need to bring in a third time point to achieve unique optima.

## Understanding Markov Model Components

To understand how the time and state dependency parameter choices influence the SOPs, let's fix $\alpha$ at the values above and vary the other parameter types one at a time.  Define a function $g$ that computes the linear predictor for non-intercept terms.  The `h` function defines the default settings for parameters not currently being varied.

```{r studyparms}
times <- c(1, 3, 7, 14, 28)
g <- function(yprev, t, gap, X, parameter=-0.5, extra) {
	tau   <- extra[1:2]
	gamma <- extra[3]
	kappa <- extra[4]
	(tau[1] * (yprev == 1) + 0*tau[2] * (yprev == 2) + tau[2] * (yprev == 3)) *
		exp(- gamma * (gap - 2)) + 
	kappa * t + parameter * (X == 2) * (t - 1) / 27
}
g <- function(yprev, t, gap, X, parameter=-0.5, extra) {
	tau   <- extra[1:2]
	gamma <- extra[3]
	kappa <- extra[4]
	tau[1] * (yprev == 1) + 0*tau[2] * (yprev == 2) + tau[2] * (yprev == 3) *
		exp(- gamma * (gap - 2)) + 
	kappa * t + parameter * (X == 2) * (t - 1) / 27
}

alpha <- qlogis(c(0.95, 0.25, 0.01))

h <- function(tau=c(-2, 0, 2), gamma=0.1, kappa=0.15, i=1) {
	extra <- c(tau, gamma, kappa)
	z <- soprobMarkovOrd(1:4, times, initial=2, absorb=4,
					intercepts=alpha, g=g, X=1, extra=extra)
	# Convert from matrix to tall and thin data frame
	data.frame(p=as.vector(z), t=as.numeric(rownames(s)[as.vector(row(s))]),
		y=as.vector(col(s)), i=i)
}

gr <- function()
	ggplot(w, aes(x=factor(t), y=p, fill=factor(y))) +
				geom_col() + facet_wrap(~ label) + 
				guides(fill=guide_legend(title='Y')) +
				xlab('t') + ylab('Occupancy Probability')

gams <- c(0.02, 0.2, 0.5, 1)
labs <- paste0('gamma=', gams)
w <- NULL
for(i in 1 : length(gams)) w <- rbind(w, h(gamma=gams[i], i=i))
w$label <- factor(w$i, labels=labs)
gr()

kappas <- c(0, .1, .2, .5)
labs <- paste0('kappa=', kappas)
w <- NULL
for(i in 1 : length(kappas)) w <- rbind(w, h(kappa=kappas[i], i=i))
w$label <- factor(w$i, labels=labs)
gr()

taus <- list(c(-2, 0, 2), c(-2, 0, 0), c(0, 0, 2), c(-4, 0, 4))
labs <- paste0('tau=', sapply(taus, function(x) paste(x, collapse=',')))
w <- NULL
for(i in 1 : length(taus)) w <- rbind(w, h(tau=taus[[i]], i=i))
w$label <- factor(w$i, labels=labs)
gr()
```

Some interpretations from these three graphs are

* Varying $\gamma$: Larger values indicate less relevance of previous state.  These don't matter very much for days 1, 3, 7 with gams 2 and 4, but for day 14 and 28, less dependence on previous state puts more patients in the highest state.
* Varying $\kappa$: The general time trend is the easiest to understand.  When $\kappa=0$ the state occupancy probabilities are stable over time except for increasing probability in the absorbing state Y=4.
* Varying $\tau$: We have set $\gamma=0.1$ so the discounting factor for the increasing measurement times with their increasing gaps is `r round(exp(-0.1 * (c(2, 4, 7, 14) - 1)), 2)`.  The $\tau$ vector captures the influence of previous state when the discounting factor is not very small, so here previous state is about $\frac{3}{4}$ ignored at day 28.  Examining the top two panels, the effect of removing the impact of previous state $y=3$ resulted in a lower probability of being in state 4.  For the bottom left panel, the effect of restoring weight to previous state of 3 and removing weight from previous state 1 resulted in little change from the first $\tau$ configuration.  In the bottom right panel, giving more weight to previous states of 1 and 3 resulted in a higher change of Y=1 on days 3 and 7, a lower chance of being at Y=2 on day 14, and a lower chance of being at Y=2 or 3 on day 28, with an appreciable increase in the chance of Y=4.

Keep in mind that these interpretations are in the context of the other parameters not currently being varied, i.e., the general time trend, intercepts, and decay function in the last graph.

## Specification of Other Parameters and Use of `Hmisc` Functions

The `Hmisc` functions we will use require a user-specified function `g` like the one above that computes the transition model linear predictor other than the proportional odds model intercept terms.  `g` computes the linear predictor for one observation, or for a series of initial states (Y values at baseline = time 0).  The arguments to `g` are the previous value of Y (the initial state if `t=1`), the current time `t`, the gap between the previous measurement time and `t`, covariate setting `X`, and `parameter` which specifies the true treatment effect.  Other parameters such as the dependence structural parameters are passed as a vector `extra`.

The transition model was defined as `g` above.  Let's use the intercept values $\alpha$ as starting estimates in the context of solving for the entire set of parameters to meet our criteria.  The `Hmisc` `intMarkovOrd` function uses the standard R optimization function `nlm` to compute the intercepts and other parameters satisfying given occupancy probability targets, once the user specifies initial guesses for all parameters.  Note that the intercept values must be in decreasing order.  `nlm` uses an efficient trial-and-error process to compute the parameter values that provide the best compromise solution.  The criterion being minimized is the sum of absolute differences between all target probabilities and the actual achieved probabilities.  With measurements made on days 1, 3, 7, 14, 28, our target values are state occupancy probabilities of 0.7, 0.14, 0.1, and 0.05 for, respectively, Y=1, 2, 3, 4 at day 28, as used previously.  These probabilities pertain to patients who start with Y=2 at baseline in treatment group 1.  Recall that he weighting function for time gaps is taken to be a negative exponential.  The decay constant is estimated during the optimization.  We apply a constraint during the optimization, namely that the decay constant is positive.

```{r compints}
target <- rbind(
	'1'  = c(0.05, 0.70, 0.24, 0.01),
	'28' = c(0.70, 0.18, 0.07, 0.05) )
extra <- c(-2, 2, 0.1, -0.15)
#extra <- c(-2, .16, 2.1, 0, 0)
# Try different starting values to see if the optimizer finds different solutions
#extra <- c(-2, 0, 2, -.15, 0)
con <- function(ints, extra) extra[3] >= 0.
for(se in 15:20) {
	cat('seed:', se, '\n')
set.seed(se)
ntry <- 100
m <- numeric(100)
ex <- matrix(0, nrow=ntry, ncol=4)
for(i in 1 : ntry) {
	u <- runif(length(extra), -1, 1)
	ex[i, ] <- extra + u
	m[i] <- intMarkovOrd(1:4, times, initial=2, absorb=4,
                  intercepts=alpha, g=g, target=target, t=c(1, 28), extra=extra + u,
									onlycrit=TRUE, constraints=con)
}
cat('Minimum sum of absolute errors:', min(m), '\n')
extra <- ex[which.min(m), ]
cat('Extra achieving this minimum:', extra, '\n')
}

z <- intMarkovOrd(1:4, times, initial=2, absorb=4,
                  intercepts=alpha, g=g, target=target, t=c(1, 28), extra=extra,
									constraints=con)
# 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(g)$extra <- z$extra
# Save intercepts
ints <- z$intercepts

# Reproduce the result using soprobMarkovOrd directly
s <- soprobMarkovOrd(1:4, times, initial=2, absorb=4,
                     intercepts=ints, g=g, X=1)  # extra already stored in function, extra=extra)
round(s, 3)

# Plot the computed general time trend
kappa  <- extra[4:5]
tt <- kappa[1] * times + kappa[2] * (times - 14) ^ 2
plot(times, tt, type='b', xlab='t', ylab='Time Effect')
```

The optimizer could achieve the desired contraints over the first and last time points.  The estimate of $\gamma$ is in the wrong direction (though is small), so $gamma$ is constrained to be zero by removing it from the parameter list.

```{r nogamma}
g <- function(yprev, t, gap, X, parameter=-0.5, extra) {
	tau   <- extra[1:3]
	kappa <- extra[4:5]
	tau[1] * (yprev == 1) + tau[2] * (yprev == 2) + tau[3] * (yprev == 3) +
	kappa[1] * t + kappa[2] * (t -14) ^ 2 + parameter * (X == 2) * (t - 1) / 27
}
extra <- c(-2, 0, 2, -0.15, 0)
alpha2 <- c(3, -1, -4.5); alpha2 <- alpha
extra <- c(-3, 0, 3, -.15, 0)
z <- intMarkovOrd(1:4, times, initial=2, absorb=4,
                  intercepts=alpha2, g=g, target=target, t=c(1, 28), extra=extra)
# 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(g)$extra <- z$extra
# Save intercepts
ints <- z$intercepts


# Repeat without absorbing state
aextra <- c(-3, 0, 3, .08, .05)
zna <- intMarkovOrd(1:4, times, initial=2,
                    intercepts=c(3, -1, -3.5), g=g, target=target, t=c(1, 28), extra=aextra)
```

Make sure the sum of absolute errors is small, otherwise the optimization algorithm may have been fooled and you'll need to try different starting values in `start`.

Before running the clinical trial simulation, simulate a large sample size for one trial to check that the simulation is working correctly.  Unlike what follows later, we will carry the absorbing state forward so that we can compute state occupancy proportions easily.

```{r checksim}
s <- simMarkovOrd(n=40000, 1:4, times, initial=2, X=c(group=1),
									absorb=4, intercepts=ints, g=g, carry=TRUE)
relfreq <- function(x) table(x) / length(x)
w       <- with(s, tapply(y, time, relfreq))
do.call(rbind, w)
```

Check the correlation structure after reshaping the dataset to have different times in columns for a subject.

```{r chkcorr}
setDT(s, key=c('id', 'time'))
w <- dcast(s, id ~ time, value.var='y')
w <- as.matrix(w)[, -1]
round(cor(w), 2)
```

# Simulation

Using the above model we simulate, for each treatment effect, 1000 clinical trials each with 600 observations.  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.  We use a distribution of initial states 1, 2, 3 and sample from that distribution to get each patient's baseline state.  Because the contrast of interest needs to take into account a treatment $\times$ time interaction, the `groupContrast` argument must be specified below.  A second set of simulations is also run and stored for later, just for OR = 0.6 and using a logistic model that has two unnecessary parameters: a square term for the time effect and a square term for the time $\times$ treatment interaction.  This second simulation will be used to see how much power is lost by allowing for more complexity involving treatment.

```{r sim}
initial <- c('1'=0.02, '2'=0.75, '3'=0.23)   # probabilities of being in baseline states
ors     <- c(0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.25)
formula  <- y ~ yprev * gap + time * group
formula2 <- y ~ yprev * gap + pol(time, 2) * group
nsim <- 1000
if(simdone) {
	est  <- readRDS('sim.rds')
	est2 <- readRDS('sim2.rds')
	} else {
	set.seed(1)
		ors <- 1.
	# Use g and ints from newer version 
	U <- readRDS('../g.rds')
	ints <- U$ints; g <- U$g
	
	# 1000*8 simulations (8 parameter values) takes about 8 minutes in all
  s <- estSeqMarkovOrd(1:4, times, initial, absorb=4, intercepts=ints,
    parameter=log(ors), looks=600, g=g, formula=formula,
    groupContrast=list(list(group=1, time=28), list(group=2, time=28)),
    timecriterion=function(y) y == 1, coxzph=TRUE, nsim=nsim, progress=FALSE)
  saveRDS(est, 'sim.rds')
  set.seed(2)
  est2 <- estSeqMarkovOrd(1:4, times, initial, absorb=4, intercepts=ints,
    parameter=log(0.6), looks=600, g=g, formula=formula2,
    groupContrast=list(list(group=1, time=28), list(group=2, time=28)),
    nsim=nsim, progress=FALSE)
  saveRDS(est2, 'sim2.rds')
}
setDT(est, key='parameter')
est[, OR := exp(parameter)]
setDT(est2)
est2[, OR := exp(parameter)]
```

Before doing the main analyses on transition model parameters, look at the computed event times and estimate the type I assertion probability $\alpha$ and power of the "time to Y=1" variable analyzed with a Cox model, and look at the the magnitude of non-proportional hazards.  Also compute the average log hazard ratio to see how the reciprocal of its antilog relates to the odds ratio (reciprocal because the Cox model is fit on time until a _good_ outcome).

```{r cox,w=6,h=3.5}
s  <- est
et <- attr(s, 'etimefreq')
et <- apply(et, 2:4, sum)
pr('Frequency distribution of times to Y=1 by group')
for(pm in dimnames(et)[[1]]) {
	cat('\nOR:', exp(as.numeric(pm)), '\n')
	print(et[pm, ,])
}
rn <- function(x) round(x, 2)
pr('Comparison of OR and HR', rn(s[, .(HR=exp(- mean(loghr))), by=OR]))
z <- rbind(data.frame(x=s[, lrchisq], Test='Cox model group comparison'),
					 data.frame(x=s[, phchisq], Test='Cox Model PH test'))
ggplot(z, aes(x=x)) + geom_histogram(bins=50) + facet_wrap(~ Test) +
	xlab(expression(chi^2))
pr('Power of Cox test for time until Y=1', 
  rn(s[, .(power=mean(lrchisq > 3.84)), by=OR]) )
pr('Proportion of simulations with chi-sq for PH > 3.84',
  rn(s[, .(rejectPH=mean(phchisq > 3.84)), by=OR]) )
```

The hazard ratio, after taking the reciprocal so that a large HR means a worse outcome, is not equatable to the odds ratio except at the null value of 1.0.

There is little evidence of non-proportional hazards.

Compute the power of the Markov model treatment comparison, and compare the $\chi^2$ statistic from the transition model to that from the Cox model.  For plotting, we take the square root of the $\chi^2$ statistics to get a more symmetric distribution over the simulations.  There is one panel per odds ratio.  For OR=0.6 also show the power of a transition model that overfitted time as quadratic.

```{r coxcompare}
s[, tchi := est * est / vest]
pr('Power of Markov model test at 28d',
   s[, .(power=mean(tchi > 3.84)), by=OR] )
pr('Power of quadratic in time Markov model at 28d when OR=0.6',
   est2[, .(power=mean(est * est / vest > 3.84))] )
rho <- s[, spearman(tchi, lrchisq)]
pr('Spearman rho correlation between Markov and Cox model chi-square:',
   inline=rn(rho))
gr <- function(x, y, z) 
	ggfreqScatter(x, y, ylab='PO Model Z Statistic', xlab='Cox Model Z Statistic',
								by=z, bins=25) +
		  geom_abline(intercept=0, slope=1, col='red', alpha=0.3) +
	    theme(legend.position='bottom', axis.text.x = element_text(angle = 90))
s[, gr(sqrt(lrchisq), sqrt(tchi), exp(parameter))]
```

The power of the longitudinal transition model is superior to that of the time to Y=1 comparison.
The $\chi^2$ statistics from the Markov ordinal longitudinal model are seen to dominate those from the Cox model for time until Y=1 as judged by the scatterplot.  The Spearman $\rho$ between the two statistics is `r rn(rho)`.  There is some loss of power to detect OR=0.6 when the quadratic model is substituted for the true linear (in time) model.


# Bayesian Power

Assume that the distribution of the maximum likelihood estimates is approximately Gaussian, and use a Gaussian prior distribution for the parameters of interest.  Then we can use the maximum likelihood estimates already simulated to get approximate Gaussian Bayesian posterior distributions, and quickly compute such things as Bayesian power, e.g., the probability that the posterior probability of a beneficial effect exceeds 0.95 at the end of the study.

Define the Bayesian assertions and priors to be used for them

* Assertion 1: log(OR) < 0 under prior with prior mean 0 and sigma: P(OR>2)=0.025
* Assertion 2: log(OR) < 0 under flat prior
* Assertion 3: log(OR) > 0 under flat prior (sigma=100)
* Assertion 4: log(OR) > 0 under optimistic prior with mean log(0.85), sigma=0.5

```{r bayespost}
asserts <- list(list('Efficacy',                   '<', 0, cutprior=log(2), tailprob=0.025),
                list('Efficacy flat',              '<', 0, mu=0,         sigma=100),
                list('Inefficacy/harm flat',       '>', 0, mu=0,         sigma=100),
                list('Inefficacy/harm optimistic', '>', 0, mu=log(0.85), sigma=0.5))
  
s <- gbayesSeqSim(est, asserts=asserts)

head(s)
attr(s, 'asserts')
alabels <- attr(s, 'alabels')  # named vector to map p1 p2 p3 p4 to labels
```

First let's examine the effect of the priors by making two pairwise comparisons: differences in posterior probabilities of efficacy under skeptical vs. flat prior, and differences in posterior probabilities of inefficacy under flat and optimistic priors.

```{r priorimpact}
w <- data.table(s)
u <- w[, .(p12max=max(abs(p1 - p2)), p12mean=mean(abs(p1 - p2)),
           p34max=max(abs(p3 - p4)), p34mean=mean(abs(p3 - p4))), by=.(look)]
z  <- melt(u, measure.vars=c('p12max', 'p12mean', 'p34max', 'p34mean'),
           variable.name='which', value.name='diff')
k <- c(p12max='Efficacy max',   p12mean='Efficacy mean', 
       p34max='Inefficacy max', p34mean='Inefficacy mean')
z[, w := k[which]]
z
```

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

```{r hits}
# Reshape results into taller and thinner data table so can plot over 3 assertions

ps <- names(alabels)
m  <- melt(w,
           measure.vars=list(ps, paste0('mean', 1:4), paste0('sd', 1:4)),
           variable.name='assert', value.name=c('p', 'mean', 'sd'))
m[, assert := alabels[assert]]
head(m)

# Define targets
target <- c(Efficacy                     = 0.95,
            'Efficacy flat'              = 0.95,
            'Inefficacy/harm flat'       = 0.9,
            'Inefficacy/harm optimistic' = 0.9)

m[, target := target[assert]]	  # spreads targets to all rows
# hit = 0/1 indicator if hitting target at the single fixed sample size
u <- m[, .(hit = mean(p > target)), by=.(OR, assert)]
u
ggplot(u, aes(x=OR, y=hit, color=assert)) + geom_line() +
  xlab('OR') + ylab('Proportion Hitting Posterior Probability Target') +
  guides(color=guide_legend(title='Assertion')) 
```

# Frequentist Power of Single Day Ordinal Outcomes

Let's consider only the detection of OR=0.6 for an unadjusted proportional odds two-group comparison of ordinal outcomes measured at a single day.  When simulating the data using the same model as above, we carry absorbing states forward here.  So once Y=4 occurs on a given day, Y=4 is considered to be in effect at all later days.

```{r single}
nsim <- 1000
if(simdone.single) ch <- readRDS('simch.rds') else {
	ch <- matrix(NA, nrow=nsim, ncol=length(times))
	colnames(ch) <- as.character(times)
	set.seed(3)
	for(isim in 1 : nsim) {
		s1 <- simMarkovOrd(n=300, 1:4, times, initial=2, X=c(group=1),
											 absorb=4, intercepts=ints, g=g, carry=TRUE)
		s2 <- simMarkovOrd(n=300, 1:4, times, initial=2, X=c(group=2),
											 absorb=4, intercepts=ints, g=g, carry=TRUE)
		s <- rbind(s1, s2)
		for(tim in times) {
			f <- lrm(y ~ group, data=s, subset=time == tim)
			ch[isim, as.character(tim)] <- f$stats['Model L.R.']
		}
	}
	saveRDS(ch, 'simch.rds')
}
pr('Frequentist power of individual day comparisons',
   apply(ch, 2, function(x) mean(x > 3.84)) )
```

The power for testing differences on day 1 has to only be $\alpha$ because the true treatment effect is zero on that day.  The power increases as time marches on.  But even on day 28 the power is significantly below the power of the ordinal longitudinal model that uses all days.

# More Information

* Full R markdown script: [hbiostat.org/R/Hmisc/simMarkovOrd.Rmd](https://hbiostat.org/R/Hmisc/simMarkovOrd.Rmd)
* COVID-19 statistical resources: [hbiostat.org/proj/covid19](https://hbiostat.org/proj/covid19)
* Bayesian design and analysis resources: [hbiostat.org/bayes](https://hbiostat.org/bayes)
* [Source code](https://github.com/harrelfe/Hmisc/blob/master/R/simMarkovOrd.r) for relevant `Hmisc` functions

# Computing Environment

`r markupSpecs$html$session()`