State:', y, '

t=', t, '

p=', round(p, 3))) # u$txt <- with(u, paste0('Group:', X, '

', label, # '

Y(3):', yp, # '

Y(7):', y, '

p=', round(p, 3))) # Graph would not simultaneously handle color and linetype # with ggplotly. Use regular ggplot object. stayc <- paste0('. P(Y(14)=1 | Y(7)=1) = ', paste(round(stay1, 2), collapse=', '), '.') caps <- paste(c('State transition probabilities', 'State occupancy probabilities'), cap, stayc) g1 <- ggplot(u, aes(x=yp, y=p, color=factor(y), linetype=factor(X))) + geom_line() + facet_wrap(~ label) + guides(color = guide_legend(title='Y'), linetype = guide_legend(title='Group')) + xlab('State on Day 3') + ylab('Probability on Day 7') g2 <- ggplot(w, aes(x=factor(t), y=p, fill=factor(y), label=txt)) + geom_col() + facet_wrap(~ label) + guides(fill=guide_legend(title='Y')) + xlab('t') + ylab('Occupancy Probability') r <- list(g1, ggp(g2)) # Create unique chunk names even if vary1 called multiple times # If omit chunk names, LaTeX -> pdf screwed up all the figures cnames <- paste0('studyparmsvary1', letters[(nvar1 + 1) : (nvar1 + 2)]) nvar1 <<- nvar1 + 2 markupSpecs$html$mdchunk(robj=r, caption=caps, cnames=cnames) invisible() } nvar1 <<- 0 cap <- 'varying $\\\\gamma$, with $\\\\tau=-4, 2$, $\\\\kappa=-0.1, 0, 0$' gams <- list(c(-1, 0), c(-1, 1), c(0, 1), c(1, 1)) vary1(gamma=gams, cap=cap) cap <- 'varying $\\\\kappa$, with $\\\\tau=-4, 2$, $\\\\gamma=1, -1$' kappas <- list(c(0, 0, 0), c(0.15, 0, 0), c(-.15, 0, 0), c(.15, .1, 0), c(.15, .1, .1)) vary1(kappa=kappas, cap=cap) cap <- 'varying $\\\\tau$, with $\\\\gamma=1, -1$, $\\\\kappa=-0.1, 0, 0$' taus <- list(c(-2, 2), c(-2, 0), c(0, 2), c(-4, 4)) vary1(tau=taus, cap=cap) ``` Keep in mind that as defined in the `h` function above, the parameter values used for the parameters not being varied in a given graph must be taken into account in interpreting the graphs. The effects of one type of parameter are most easily seen in the transition probability graphs. ## 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 `yprev` 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. `g` returns a matrix with number of rows equal to the length of `yprev` and number of columns equal to one if the model is PO, or equal to the number of intercepts if it is partial PO. Other parameters such as the dependence structural parameters are passed as a vector `extra`. When non-proportional odds is modeled, as in the example above, `g` returns more than one column corresponding to all the intercept increments coming from the partial proportional odds model. The transition model was defined as `g` above. Let's use the previous 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 second set of target values are 28d state occupancy probabilities of 0.7, 0.14, 0.1, and 0.05 for, respectively, Y=1, 2, 3, 4. 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. Define a function `ftarget` to specify the $y=1$ transition probability constraint. ```{r compints,fig.cap='State occupancy probabilities with initial state 2 and absorbing state 4. Each pair of bars represents treatment 1 (left) and treatment 2 (right).'} # State occupancy probability targets target <- rbind( # row names define the times for which constraints apply '1' = c(0.05, 0.70, 0.24, 0.01), '28' = c(0.70, 0.18, 0.07, 0.05) ) # Transition probability target ftarget <- function(intercepts, extra) { # Constrain P(Y(14) = 1 | Y(7) = 1, X=1) = 0.9 # Pr(Y = 1) = 1 - Pr(Y > 1) which uses first intercept p1 <- 1. - plogis(intercepts[1] + g(1, 14, 7, 1, extra=extra)[1]) abs(p1 - 0.9) } # Try different starting values z <- findstart(g=g, extra=c(tau=c(-3, 2), gamma=c(0,0), kappa=c(kappa1, 0, 0)), intercepts=alpha, times=times, ftarget=ftarget) # Graph and optionally print (if pdf output) sop12(y=1:4, times=times, initial=2, absorb=4, intercepts=z$intercepts, g=g, extra=z$extra) ``` ```{r compints2} # 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 # Save g object and intercepts saveRDS(list(ints=ints, g=g), 'g.rds') ``` Let's check our $Y(7)$ to $Y(14)$ transition probability to see to what extent it is satisfied. ```{r check1} ftarget(ints, z$extra) xb <- ints[1] + g(1, 14, 7, 1)[1] 1. - plogis(xb) ``` ```{r repro,fig.cap='State occupancy probabilities with initial state 2 and absorbing state 4, for group 1'} # 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 plotsop(s) ``` Since our model contains departures from proportional odds, check that the implied exceedance probabilities are in descending order as $y \uparrow$. Only time interacts with the intercept increments, so we only need to vary time. ```{r checkprobs} chkints(times=times, intercepts=ints, extra=formals(g)$extra) ``` We also estimate simulation model parameters when there is no absorbing state, just for comparison. ```{r noabs,fig.cap='State occupancy probabilities for initial state 2 without an absorbing state'} zna <- intMarkovOrd(1:4, times, initial=2, intercepts=c(4, 0, -3), g=g, target=target, t=c(1, 28), ftarget=ftarget, extra=c(-1, 3, .5, 2, -.5, 0, -.5)) sop12(y=1:4, times=times, initial=2, intercepts=zna$intercepts, g=g, extra=z$extra) ``` 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 for `intercepts` and/or `extra`. 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. A graphic shows the within-patient correlation matrix of ordinal outcomes across time. ```{r checksim,fig.cap='Correlations between numeric states at four measurement times within 10,000 simulated patients with initial state 2 and absorbing state 4'} r <- testsamp(intercepts=ints, times=times) w <- plotCorrM(r) ggp(w[[1]]) ``` ```{r checksim2,fig.cap='Correlation vs. time gap, computed from 10,000 simulated patients'} ggp(w[[2]]) ``` The correlation pattern from the Markov simulated data is a typical serial correlation pattern with waning correlation as the time gap expands. Let's also simulate a single 10000-patient trial to compare the usual proportional odds logistic model parameter variances with those from the cluster sandwich robust estimator. ```{r simse} n <- 10000 set.seed(8) s1 <- simMarkovOrd(n=n / 2, y=1:4, times, initial=2, X=c(group=1), absorb=4, intercepts=ints, g=g) s2 <- simMarkovOrd(n=n / 2, y=1:4, times, initial=2, X=c(group=2), absorb=4, intercepts=ints, g=g) s2$id <- s2$id + n s <- rbind(s1, s2) s$yprev <- as.factor(s$yprev) s$group <- as.factor(s$group) dd <- datadist(s); options(datadist='dd') f <- lrm(y ~ yprev * pmax(gap - 2, 0) + time * group, data=s, x=TRUE, y=TRUE) frobust <- robcov(f, s$id) # Define needed group contrast at 28 times since time x group interaction present g1 <- list(group=1, time=28); g2 <- list(group=2, time=28) contrast(f, g2, g1) contrast(frobust, g2, g1) ``` The estimated standard error of the 28d group effect is larger with the cluster sandwich estimate. ## Simulation Models vs. Analysis Models The statistical model used on real data is not privy to all the assumptions made by the simulation model. So typically the analysis model needs to have more parameters to account for what we don't know. Here are two examples: * In our simulation model we don't have a main effect for treatment as we assume a zero treatment effect at the first follow-up time $t=1$. The analysis model needs the treatment main effect to be added to such a treatment $\times$ time interaction effect. * Our simulation model assumes that the impact of time gap is absent when the previous state is $y=1$. But the fitted model allows for a main effect of $\max(\mathrm{gap} -2, 0)$, this main effect applying to the reference cell of previous $y=1$. Let's compare the parameter values used in the simulation with the estimated values from a partial proportional odds model fit, using the 10,000 patient simulated dataset derived just above. Take into account that the simulation model has a one-day offset for $t$. ```{r chksimmod} g # show simulation model ints s$tim <- s$time - 1 formula <- y ~ yprev * pmax(gap - 2, 0) + tim * group f <- VGAM::vgam(formula, VGAM::cumulative(parallel = FALSE ~ tim, reverse=TRUE), data=s) k <- coef(f) k # vgam parameterizes partial proportional odds terms as total effects and not as offsets # from the main proportional odds term. Rephrase to compare with our notation. kn <- c('time:1', 'time:2', 'time:3') kn <- c('tim:1', 'tim:2', 'tim:3') kappa <- k[kn] k[kn] <- c(kappa[1], kappa[2] - kappa[1], kappa[3] - kappa[1]) # From simulation model: p <- formals(g)$extra model <- c(ints, p[c('tau1', 'tau2')], 0, p[c('kappa1','kappa2','kappa3')], 0, p[c('gamma1','gamma2')], -0.5/27) compp <- round(cbind(Estimated=k, True=model), 3) compp ``` # 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 PO 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. A linear gap decay function is used. It would have perhaps been more reasonable to use an exponential decay such as $\exp(-\lambda \delta)$ but we cannot fit models that are nonlinear in the unknown parameters using standard regression software. **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` function is used to fit the partial proportional odds model. When doing Bayesian analysis, the `rmsb` package `blrm` function will fit a constrained partial PO model. When a PO model was fitted but non-PO was present, one simulation with OR=1.0 (not shown here) resulted in $\alpha > 0.1$. On an 8-core machine we use 7 cores in parallel to speed up simulations by almost a factor of 7. ## Partial Proportional Odds Model Simulation ```{r sim} initial <- c('1'=0.02, '2'=0.75, '3'=0.23) # probabilities of being in baseline states formula <- y ~ yprev * pmax(gap - 2, 0) + time * group ppo <- ~ time # partial PO model with PO exception for time ors <- c(0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.25) contrvgam <- function(fit) { # for customized vgam contrasts if non-PO in effect beta <- coef(fit) v <- vcov(fit) # Treatment effect at 28 days (linear time x treatment interaction) me <- 'group2'; ia <- 'time:group2' # names of main effect and interaction parameters effect <- beta[me] + 28 * beta[ia] if(ia %nin% rownames(v)) saveRDS(fit, 'problemfit.rds') vr <- v[me, me] + 28. * 28. * v[ia, ia] + 2. * 28. * v[me, ia] list(Contrast=effect, SE=sqrt(vr)) } sims <- dosim(times=times, g=g, initial=initial, intercepts=ints, ors=ors, formula=formula, ppo=ppo, groupContrast=contrvgam, timecriterion=function(y) y == 1, seed=4, file='simppo.rds') # Don't run this until figure out how to write a contrvgam function for it # formula2 <- y ~ yprev * pmax(gap - 2, 0) + pol(time, 2) * group # est2 <- dosim(file='sim2.rds', g=g, initial=initial, intercepts=ints, # ors=ors, formula=formula2, ppo=ppo, # groupContrast=contrvgam, timecriterion=function(y) y == 1, # seed=5) ``` 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 anti-log relates to the odds ratio (reciprocal because the Cox model is fit on time until a _good_ outcome). Then 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. ```{r analyzesim, results='asis'} examSim(sims, desc='Partial PO model with absorbing state 4') ``` 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. 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. ## Proportional Odds Model Simulation The data were simulated using a realistic partial PO model to relax how the mix of events can change over time. What happens when we improperly fit a model that excludes departures from PO with respect to time? We repeat the above simulations to find out. We use the cluster sandwich robust covariance estimator to make up for some of the lack of fit. ```{r polof,results='asis'} # Define a contrast that works with rms::lrm so that we get treatment # effect at day 28 contr <- list(list(group="1", time=28), list(group="2", time=28)) po <- dosim(times=times, g=g, initial=initial, intercepts=ints, ors=ors, formula=formula, groupContrast=contr, cscov=TRUE, seed=4, file='simpo.rds') ``` ```{r polof2,results='asis'} examSim(po, FALSE, FALSE, FALSE, FALSE, desc='PO model with lack of fit') ``` Without using a robust covariance estimate, type I $\alpha$ assertion probability was 0.14 (not shown). The robust covariance estimate fixed the $\alpha$ problem, and the group comparison assuming PO had the same power as without assuming PO. The problems with the after-the-fit variance corrections are (1) we can't be certain that the magnitude of the group effect is correct, and more obviously (2) the approach does not transport to our Bayesian model. ```{r quad, echo=FALSE, eval=FALSE} pr('Power of quadratic in time Markov model at 28d', est2[, .(power=mean(est * est / vest > 3.8415)), by=OR] ) ``` # 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(sims$sim, 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 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 easily 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(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 ptarget <- c(Efficacy = 0.95, 'Efficacy flat' = 0.95, 'Inefficacy/harm flat' = 0.9, 'Inefficacy/harm optimistic' = 0.9) m[, ptarget := ptarget[assert]] # spreads targets to all rows # hit = 0/1 indicator if hitting target at the single fixed sample size u <- m[, .(hit = mean(p > ptarget)), by=.(OR, assert)] u u$txt <- with(u, paste0('OR:', OR, '

', assert, '

Hit probability:', hit)) ggp(ggplot(u, aes(x=OR, y=hit, color=assert, label=txt)) + 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} seed <- 3 nsim <- 1000 gg <- g formals(gg)$parameter <- log(0.6) file <- 'simch.rds' # See if previous simulation had identical inputs and Hmisc code hashobj <- hashCheck(simMarkovOrd, nsim, seed, times, ints, gg, file=file) hash <- hashobj$hash ch <- hashobj$result if(! length(ch)) { # Run the simulations since something changed or file didn't exist # Define function to run simulations on one core run1 <- function(reps, showprogress, core) { ch <- matrix(NA, nrow=reps, ncol=length(times)) colnames(ch) <- as.character(times) for(isim in 1 : reps) { if(isim %% 10 == 0) showprogress(isim, reps, core) s1 <- simMarkovOrd(n=300, 1:4, times, initial=2, X=c(group=1), absorb=4, intercepts=ints, g=gg, carry=TRUE) s2 <- simMarkovOrd(n=300, 1:4, times, initial=2, X=c(group=2), absorb=4, intercepts=ints, g=gg, 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.'] } } list(ch=ch) } ch <- runParallel(run1, reps=nsim, seed=seed) attr(ch, 'hash') <- hash saveRDS(ch, file, compress='xz') } pr('Frequentist power of individual day comparisons', apply(ch, 2, function(x) mean(x > 3.8415)) ) ``` 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. # Power With More Observations Per Patient Let's repeat the main result but instead of sampling on selected days between 1 and 28, sample every day. For this setup, since the time gap is a constant 1.0 we no longer need `gap` in the model. The number of ORs simulated is reduced. ```{r studyparms28,fig.cap='State occupancy probabilities for 28 consecutive days, with absorbing state 4. Pairs of bars indicate treatment 1 (left bar) and treatment 2 (right bar).'} times <- 1 : 28 g <- function(yprev, t, gap, X, parameter=-0.5, extra) { tau <- extra[1:2] kappa <- extra[3:5] lp <- matrix(0., nrow=length(yprev), ncol=3, dimnames=list(as.character(yprev), c('2','3','4'))) # 3 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 for(yp in yprev) lp[as.character(yp), ] <- tau[1] * (yp == 1) + tau[2] * (yp == 3) + (t - 1) * (kappa[1] + c(0., kappa[2], kappa[3])) + parameter * (X == 2) * (t - 1) / 27 lp } alpha <- qlogis(c(0.95, 0.25, 0.01)) # Try different starting values z <- findstart(g=g, intercepts=alpha, extra=c(tau=c(-2, 2), kappa=c(0, 0, 0)), times=times, seed=1) sop12(y=1:4, times=times, initial=2, absorb=4, intercepts=z$intercepts, g=g, extra=z$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 ``` Reproduce the result using soprobMarkovOrd directly. ```{r sp28r,fig.cap='State occupancy probabilities with initial state 2 and absorbing state 4 for treatment group 1'} s <- soprobMarkovOrd(1:4, times, initial=2, absorb=4, intercepts=ints, g=g, X=1) plotsop(s) ``` ```{r sp28,fig.cap='Correlations between times within 10,000 simulated patients with initial state 2 and absorbing state 4'} r <- testsamp(n=10000, intercepts=ints, times=times) w <- plotCorrM(r) ggp(w[[1]]) ``` ```{r sp28b,fig.cap='Correlations vs. time gap for 10,000 simulated patients'} ggp(w[[2]]) ``` ```{r sp28no,fig.cap='State occupancy probabilities with no absorbing state'} # Repeat without absorbing state aextra <- c(tau=c(-4, 3), gamma=c(0,0), kappa=c(-0.2, 0.05, 0.15)) # initial guess x <- findstart(g=g, intercepts=c(3, -1, -4), extra=aextra, ftarget=ftarget, absorb=NULL, times=times, seed=3) sop12(y=1:4, times=times, initial=2, intercepts=x$intercepts, g=g, extra=x$extra) zna <- intMarkovOrd(1:4, times, initial=2, intercepts=x$intercepts, g=g, target=target, t=c(1, 28), ftarget=ftarget, extra=x$extra) ``` Now simulate 28d of data per patient. ```{r sim28d} initial <- c('1'=0.02, '2'=0.75, '3'=0.23) # probabilities of being in baseline states ors <- c(0.5, 0.6, 0.7, 0.8, 0.9, 1.0) formula <- y ~ yprev + time * group est <- dosim(times=times, g=g, initial=initial, intercepts=ints, ors=ors, formula=formula, ppo=ppo, groupContrast=contrvgam, timecriterion=function(y) y == 1, seed=6, file='sim28.rds') ``` As before we also check the performance of a Cox two-sample test for time to $Y=1$ as well as that of the ordinal longitudinal model. We also show Kaplan-Meier cumulative incidence estimates for being discharged to home, censoring on death (hence competing risks are not really being dealt with optimally). The cumulative incidences are computed after pooling all 1000 studies' data. ```{r exam28,results='asis'} examSim(est, desc='Partial PO model with 28 days of measurements') ``` # Operating Characteristics Under Response-Dependent Sampling In COVID-19 therapeutic trials with longitudinal outcomes, many of the trials randomize patients once they arrive to hospital. While in hospital, the ordinal scale may be assessed frequently, e.g., daily. Once the patient is discharged home, even if this is not an absorbing state (patient can return to hospital if condition worsens again), it may be difficult to do daily outcome assessments. Starting with the previous simulation of 28 days of follow-up for each patient (less for patients reaching absorbing state `Y=4`), let's not sample patients as frequently once they reach `Y=1`. Specifically, the planned assessment times starting at day $t$ that the patient first reached `Y=1` will be `t+3`, `t+6`, `t+9`, ... $\min(28, T_{a})$ where $T_{a}$ is the time until $Y=4$ if $Y$ was ever 4 in the 28 days of follow-up. We construct an `rdsample` function that accomplishes this. The function returns `NULL` if `Y=1` was not reached before day 28 or if no times were dropped after reaching `Y=1`. ```{r rdsamp} rdsamp <- function(times, y) { if(! any(y == 1)) return(NULL) # Find index of first observation reaching y=1 j <- min(which( y == 1 )) if(j == 28) return(NULL) planned.times <- c(if(j > 1) times[1 : (j - 1)], seq(times[j], max(times), by=3)) if(length(planned.times) == length(times)) return(NULL) planned.times } # Test the function list(rdsamp(1, 4), rdsamp(1:3, 2:4), rdsamp(1:28, rep(3, 28)), rdsamp(1:28, c(rep(3, 27), 1)), rdsamp(1:28, c(2, 3, 1, rep(2, 25))) ) ``` Simulate 1000 clinical trials with reduced frequency of sampling post reaching `Y=1`. ```{r simrds} initial <- c('1'=0.02, '2'=0.75, '3'=0.23) # probabilities of being in baseline states ors <- c(0.7, 0.8, 1.0) formula <- y ~ yprev * pmax(gap - 2, 0) + time * group est <- dosim(times=times, g=g, initial=initial, intercepts=ints, ors=ors, formula=formula, ppo=ppo, groupContrast=contrvgam, timecriterion=function(y) y == 1, rdsample=rdsamp, seed=2, file='simrds.rds') ``` As before we also check the performance of a Cox two-sample test for time to $Y=1$ as well as that of the ordinal longitudinal model. Let's see how much power is lost by the decreased post `Y=1` sampling frequency. Here we look at only three odds ratios. The data model must once again incorporate gap discounting by interacting `gap` with previous states. ```{r examrds,results='asis'} examSim(est, timedist=FALSE, hist=FALSE, scatter=FALSE, desc='Partial PO model with response-dependent measurement times') ``` There is some power loss with OR=0.8. The loss would be expected to be smaller had $Y=1$ been more of an absorbing state, i.e, the probability of staying in that state been greater than the 0.9 used in our simulation model. # Creating a Simulation Model From a Previous Study VIOLET 2 was a randomized clinical trial of seriously ill adults in ICUs to study the effectiveness of vitamin D vs. placebo. Vitamin D was not effective. VIOLET 2 is a rich dataset having daily ordinal outcomes assessed for 28 days. Due to an apparent anomaly in ventilator use on day 28 we restrict attention to days 1-27. The original [NEJM paper](https://www.nejm.org/doi/full/10.1056/NEJMoa1911124) focused on 1078 patients with confirmed baseline vitamin D deficiency. The analyses presented here use 1352 of the original 1360 randomized patients. We fit a partial proportional odds Markov transition model, similar to the one used above, to VIOLET 2. We start with a side analysis that examines the relative efficiency of using only some of the 27 days of observation. The VIOLET 2 data are not yet publicly available. ```{r violet2} require(VGAM) dcf <- readRDS('dcf.rds') setDT(dcf, key=c('id', 'day')) # includes records with death carried forward dcf <- dcf[! is.na(bstatus) & day < 28, .( id, day, bstatus, status, treatment, ddeath)] # Since baseline status did not distinguish ARDS from on ventilator, # we must pool follow-up status likewise s <- dcf[, bstatus] levels(s) <- list('Vent/ARDS'='ARDS', 'Vent/ARDS'='On Vent', 'In Hospital/Facility'='In Hospital') dcf[, table(bstatus, s)] dcf[, bstatus := s] d <- dcf # Get day 14 status carrying death forward status14 <- d[day == 14, status, by=.(treatment, id, bstatus)] # Before truncating records at death, get the observed frequency distribution relf <- function(y) as.list(table(y) / length(y)) w1 <- d[, relf(status), by=day] setnames(w1, 'day', 'time') w1[, time := as.integer(time)] w1[, type := 'VIOLET 2 Results'] # Compute correlation matrix w <- dcast(d[, .(id, day, y=as.numeric(status))], id ~ day, value.var='y') rviolet.obs <- cor(as.matrix(w[, -1]), use='pairwise.complete.obs') # Don't carry death forward when fitting Markov models: d <- d[day <= ddeath] setkey(d, id, day) # Get distribution of ventilator/ARDS-free days observed in VIOLET 2 # Pool treatment groups vfd <- function(time, y) # time is not used as.integer(ifelse(any(y == 'Dead'), -1, sum(y != 'Vent/ARDS'))) vvfd <- d[, .(vf = vfd(day, status)), by=id] vvfd <- relfreq(vvfd[, vf]) ``` ## Estimating Efficiency From the Study The model is simplified to assume a constant treatment effect over time. We vary the number of days of observation used, and compute the variance of the estimated treatment effect for each model. Days are sampled in such a way as to approximate equal gap times between measurements so that gaps can be ignored. Efficiency is the minimum variance (variance using measurements on all 27 days) divided by the variance using a given number of almost equally-spaced days. The non-PO part of the model has 6 parameters: 2 time components $\times$ 3 states corresponding to PO model intercepts. Two time components were required to accommodate the very low number of patients who were sent home on the first day. The linear spline used here essentially grants an exception at $t=1$. ```{r daysamp,fig.cap='Effective sample size per subject with reference to using only the response at day 14 (top panel); relative efficiency of a treatment comparison with reference to using all 28 days (bottom panel). Both graphs are based on the variance of the log odds ratio for state transitions.'} V <- numeric(27) a <- 'treatmentVitamin D' # Make sampled days as evenly distributed as possible for(ndays in 1 : 27) { if(ndays == 1) { # First analysis uses an ordinary unconditional PO model on day 14 f <- vgam(ordered(status) ~ bstatus + treatment, cumulative(reverse=TRUE, parallel=TRUE), data=status14) } else { s <- round(seq(1, 27, length=ndays)) ds <- d[day %in% s] ds[, pstatus := shift(status), by=id] ds[is.na(pstatus), pstatus := bstatus] ds[, pstatus := pstatus[drop=TRUE]] # drop unused previous status Dead f <- if(ndays < 5) vgam(ordered(status) ~ pstatus + treatment + day, cumulative(reverse=TRUE, parallel=FALSE ~ day), data=ds) else vgam(ordered(status) ~ pstatus + treatment + lsp(day, 2), cumulative(reverse=TRUE, parallel=FALSE ~ lsp(day, 2)), data=ds) } V[ndays] <- vcov(f)[a,a] } if(outfmt == 'pdf') pr(obj=round(V, 5)) w <- rbind(data.frame(days=1:27, y=V[27]/V, type='Relative Efficiency'), data.frame(days=1:27, y=V[1]/V, type='Effective Sample Size Per Subject')) txt <- with(w, paste0(type, '

', round(y, 2), '

Day:', days)) ggp(ggplot(w, aes(x=days, y=y, label=txt)) + geom_line() + facet_wrap(~ type, scales='free_y', nrow=2) + xlab('Number of Days Sampled') + ylab('')) ``` Sampling 3 days is estimated to have the effect of doubling the number of patients compared to sampling only day 14. Sampling all 27 days has the effect of increasing the sample size almost 5-fold. Evidence for non-PO can be seen by computing Wald $z$ statistics for the last six terms in the model. ```{r vnonpo} round(coef(f), 4) coef(f)[7:12] / sqrt(diag(vcov(f)[7:12, 7:12])) ``` ## Simulation Model From the Study The last model fitted was for all 27 days. The `vgam` function's non-PO terms are parameterized to pertain to each category after the first (Dead). We write a `g` function for the state transitions, and save intercepts separately. ```{r vpo} k <- coef(f) ints <- k[1 : 3] extra <- k[- (1 : 3)][-3] g <- function(yprev, t, gap, X, parameter=0, extra) { tau <- extra[1:2] kappa <- extra[3:8] lp <- matrix(0., nrow=length(yprev), ncol=3, dimnames=list(as.character(yprev), c('Vent/ARDS', 'In Hospital/Facility', 'Home'))) tp <- pmax(t - 2, 0) # second time term, from linear spline knot at t=2 # 3 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 for(yp in yprev) lp[as.character(yp), ] <- tau[1] * (yp == 'In Hospital/Facility') + tau[2] * (yp == 'Home') + t * kappa[1:3] + tp * kappa[4:6] + parameter * (X == 2) lp } formals(g)$extra <- extra ``` Simulate $n$ patients on placebo and compare the state occupancy frequencies with those observed in the trial, pooling treatments. $n$ is set to $10 \times$ the number of patients in the trial sample we are using to make comparisons easy after we divide frequencies by 10. Just for comparison purposes we carry death forward. For the initial state we sample from the observed proportions. ```{r vcomp,fig.cap='State occupancy probabilities averaged over the observed VIOLET 2 distribution of initial states, for treatment 1'} n <- length(unique(dcf$id)) * 10 n initialProp <- d[day==1, relfreq(bstatus)] set.seed(13) initialSamp <- sample(names(initialProp), n, TRUE, prob=initialProp) # Compute state occupancy probabilities for group 1 # We need to compute a weighted average over the result for each baseline state # sop <- list() for(initial in names(initialProp)) sop[[initial]] <- soprobMarkovOrd(levels(d$status), 1:27, initial=initial, absorb='Dead', intercepts=ints, g=g, X=1) avgProp <- initialProp['Vent/ARDS'] * sop[['Vent/ARDS']] + initialProp['In Hospital/Facility'] * sop[['In Hospital/Facility']] w2 <- data.table(avgProp) w2[, time := as.numeric(rownames(w2))] w2[, type := 'Model Fitted to VIOLET 2'] plotsop(avgProp) s <- simMarkovOrd(n=n, y=levels(d$status), times=1:27, initial=initialSamp, X=c(group=1), absorb='Dead', intercepts=ints, g=g, carry=TRUE) # Relative frequencies just for day 27 with(subset(s, time == 27), relfreq(y)) with(subset(dcf, day == 27), relfreq(status)) s <- setDT(s, key=c('id', 'time')) w3 <- s[, relf(y), by=time] w3[, type := 'Data Simulated from Fitted Model'] w <- dcast(s[, .(id, time, y=as.numeric(y))], id ~ time, value.var='y') rviolet.sim <- cor(as.matrix(w[, -1])) ``` ## Compare Serial Correlation Patterns from the Study and Simulated Data Let's look at a variogram-like display of the observed correlations between times within subject to see how correlations vary with changing gap in measurement times. We compare this to the correlations seen in the simulated data. ```{r violetr,fig.cap='Serial correlation pattern within patient over time observed in VIOLET 2 along with correlations computed from a large random sample from the fitted partial proportional odds ordinal state transition model. In this variogram, the curves going horizontally and downwards represent individual days, and the vertical spread represents a non-isotonic serial correlation pattern. In other words, the correlation between two times within patient depends not only on the (unsigned) time gap but also on absolute time.'} dv <- cbind(plotCorrM(rviolet.obs, 'data'), type='Observed') ds <- cbind(plotCorrM(rviolet.sim, 'data'), type='Simulated') db <- rbind(dv, ds) ggp(ggplot(db, aes(x=delta, y=r, label=txt)) + geom_point() + geom_smooth(method=loess) + facet_wrap(~ type) + xlab('Absolute Time Difference, Days') + ylab('r')) ``` The model's correlation pattern is in good agreement with that observed in the study data, except for day 1. ## Compare State Probabilities from Study, Model Parameters, and Data Simulated from the Model ```{r compp,fig.cap='Observed VIOLET 2 state occupancy probabilities, probabilities computed exactly from the fitted model parameters, and proportions from simulating new patients from the fitted model'} w <- rbind(w1, w2, w3) r <- melt(w, measure.vars=levels(dcf$status), variable.name='y', value.name='p') txt <- with(r, paste0(type, '

', y, '

t=', time, '

p=', round(p, 3))) ggp(ggplot(r, aes(x=time, y=p, fill=y, label=txt)) + geom_col() + facet_wrap(~ type, ncol=1) + guides(fill = guide_legend(title='')) + xlab('Day') + ylab('Probability')) ``` ```{r compp2,fig.cap='Same quantities as in the previous graph but separated by time so that agreement across the three estimates may be more readily seen'} gt <- if(outfmt == 'pdf') 'Simulated data, fitted model, and VIOLET 2 results (left to right)' else '' ggp(ggplot(r, aes(x=type, y=p, color=y, label=txt)) + geom_point() + facet_wrap(~ time) + guides(color = guide_legend(title='')) + theme(axis.text.x = element_blank()) + # element_text(angle=90)) + xlab('') + ylab('Probability') + ggtitle(gt)) ``` ## Compare Distribution of Ventilator/ARDS-free Days from Study and Data Simulated from the Model Here we use data simulated under OR=1. ```{r compvfd,fig.cap='Distribution of observed VIOLET 2 ventilator/ARDS-free days (-1 for death) compared with the same summary measure computed on data simulated from the model'} dv <- data.frame(y=as.numeric(names(vvfd)), unname(vvfd), type='VIOLET 2') # Summarize simulated data simvf <- s[, .(y = vfd(time, y)), by=id] r <- relfreq(simvf$y) d2 <- data.frame(y=as.numeric(names(r)), unname(r), type='Data Simulated from Model') dv <- rbind(dv, d2) dv$txt <- with(dv, paste0(type, '

', round(Freq, 3))) ggp(ggplot(dv, aes(x=y, y=Freq, label=txt)) + geom_col() + facet_wrap(~ type, nrow=2) + guides(color = guide_legend(title='')) + xlab('Ventilator/ARDS-free Days') + ylab('Proportion')) ``` ## Simulate Performance of Statistical Methods on Studies Like VIOLET 2 Next simulate 1000 clinical trials using this VIOLET 2-based simulation model, with treatment odds ratios of 1.0 and 1.3 (it's greater than 1 because larger $Y$ are better in this setting) for ordinal state transition probabilities. The main purposes of this simulation is to compare the frequentist power of the Markov PO model with that of a "time to home" Cox model test, and with a Wilcoxon test of ventilator/ARDS-free days counting death as -1 day. The clinical trials, unlike VIOLET 2, have 600 patients. ```{r vsim} formula <- ordered(y) ~ yprev + group + lsp(time, 2) ppo <- ~ lsp(time, 2) est <- dosim(y=levels(d$status), absorb='Dead', times=1:27, g=g, initial=c(Home=0, initialProp), intercepts=ints, ors=c(1, 1.3), formula=formula, ppo=ppo, timecriterion=function(y) y == 'Home', sstat=vfd, seed=13, file='simv.rds') ``` ```{r vsim2,results='asis'} examSim(est, hist=FALSE, timeevent='home', desc='Simulation of VIOLET 2') ``` # 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()`