Ordinal Regression with Longitudinal Model
Berry Consultants, LLC
1 Introduction
We consider here a trial with an ordinal primary endpoint, and an underlying longitudinal model that uses information from previous visits for all subjects who have not reached the primary endpoint visit. For illustration purposes, we consider a trial with:
- control arm
- treatment arm
The ordinal scale ranges from \(0\) to \(5\), where \(5\) is the worst possible outcome. Suppose that the primary model analyzes the ordinal endpoint at 180 days, but that we also have that same endpoint measured at 30 and at 90 days. We will use a longitudinal model to include information from those early visits.
2 Ordinal Logistic Regression with Markov Longitudinal Model
2.1 Ordinal Regression
Ordinal data are a special case of categorical data where the ordering of the discrete categories matters. While categorical data with unordered categories are often modeled with an exchangeable probability model for the latent category probabilities, such an assumption is not desirable in ordinal data. In fact, the ordering of the categories implies a certain degree of correlation, i.e., neighboring categories should be more correlated than distant ones.
A common way to model discrete ordinal data is to use a generative model for a continuous latent variable that is then censored on the set of discrete possible values to yield the ordinal probabilities. We begin by introducing a latent continuous space such as the real line, \(X = \mathbb{R}\), with a probability distribution specified by the probability density function (p.d.f.) \(\pi(x)\) and associated cumulative distribution function (c.d.f.) \(\Pi(x)\). In the following, we will use a c.d.f. given by the logistic function, i.e. \[\Pi(x) = \frac{1}{1+e^{-x}},\] which has the associated p.d.f. \[\pi(x)= \frac{e^{-x}}{(1+e^{-x})^{2}}.\] This choice yields an ordinal logistic model. We could have also chosen a Gaussian probability density, which defines an ordinal probit model. In the figures below, we illustrate the c.d.f. and p.d.f. that define the ordinal logistic model.
Next, to model \(K\) ordered classes we partition \(X\) into \(K\) intervals using \(K+1\) cut points, \(\{c_{0} = -\infty, c_{1}, \dots, c_{K-1}, c_{K} = +\infty\}\), from which we can compute \(K\) ordinal probabilities as differences between the cumulative distribution function evaluated at the cut points, i.e. \[p_{k}= \Pi(c_{k}) - \Pi(c_{k-1}), \ k = 1, \dots, K.\] The partitioning of \(X\) ensures that the category probabilities sum to one, and the ordering of cut points induces the desired correlations between ordinal probabilities.
The two figures below illustrate how the ordinal probabilities are defined starting from the distribution \(\Pi(x)\) (or, equivalently, \(\pi(x)\)) and the set of cut points \(c_{0}, \dots, c_{K}\). Looking at the p.d.f. \(\pi(x)\), the ordinal probabilities \(p_{k}\) are defined as the area under the curve corresponding to each interval \([c_{k-1}, c_{k}]\). We will treat the cutpoints as model parameters that affect on a latent scale the response probabilities. For instance, lowering the cut point \(c_{k}\) lowers the probability allocated to the \(k^{th}\) ordinal while increasing the probability allocated to the neighboring \((k+1)^{st}\) ordinal. Because of their fixed values the bounding cut points, \(c_{0} = -\infty\) and \(c_{K} = +\infty\), are sometimes ignored in practice such that the cut points are defined as only the \(K-1\) interior boundaries.
To model observed counts in each category, we can complete the model with a Multinomial observational model \[n_{1}, \dots, n_{K} \mid (p_{1}, \dots, p_{K}) \sim \text{Multinomial} (p_{1}, \dots, p_{K}).\]
To include a treatment effect that differentiates the ordinal probabilties between arms, we introduce a latent parameter \(\theta\) that can “drag” the probability masses to higher or lower ordinals. The effect of this parameter is usually taken as a linear effect on the latent scale, i.e. by translating all the cutpoints from \(c_{j}\) to \(c_{j} + \theta\). When \(\theta\) is positive this translation is equivalent to shifting the latent probability density function to smaller values, concentrating the integrated probabilities to smaller ordinals as desired. Similarly, when \(\theta\) is negative this translation is equivalent to shifting the latent probability density function to larger values, concentrating the integrated probabilities to larger ordinals.
Consequently, in the presence of a treatment effect the ordinal probabilities become
\[p_{k} = \Pi(c_{k} + \theta) - \Pi(c_{k−1} + \theta).\]
Using the properties of the logistic function, we can rewrite the probabilities as
\[p_{k} = \Pi(- c_{k−1} - \theta) - \Pi(- c_{k} - \theta),\]
which is for example how the ordered_logistic_lpmf(0, c + theta) function is defined in the Stan functions reference guide.
For instance, \(e^{\theta} = 1.8 > 1\) has the effect of moving the cutpoints to larger values and, hence, the probability masses towards smaller values, which in our example corresponds to patient benefit (smaller values being more beneficial to patients). Below are the updated cutpoints and the resulting ordinal probabilities.
2.2 Longitudinal Model
At the time of each analysis there will be subjects who have an unknown 180-day response value. We use the 30-day and 90-day response values as possibly predictive of the 180-day response, allowing subjects with these earlier measurements to be included in the analyses of the 180-day measurement. The longitudinal model allows for learning the relationship between the 30-day and 90-day response values as the accruing empirical data is used to determine the strength of the association between the values for each treatment group. Analyses of the 180-day response values are performed with multiple imputation from the longitudinal model for patients with an unknown 180-day response value.
In the case of ordinal outcomes, a possible strategy is to assume a Markov structure for the transition between visits, i.e., assume that the response value at a certain visit only depends on the response value at the previous visit. The initial 30-day responses are assigned Dirichlet distributions. In particular, we assume a Dirichlet prior for the vector of probabilities \(\mathbf{p}^{(30)}\) corresponding to each response value at 30 days for each treatment arm separately, i.e., \[ \begin{aligned} &(p_{ctr,0}^{(30)}, \dots, p_{ctr,5}^{(30)} ) \sim \text{Dirichlet}(\beta_{0}, \dots, \beta_{5}), \\ &(p_{trt,0}^{(30)}, \dots, p_{trt,5}^{(30)} ) \sim \text{Dirichlet}(\beta_{0}, \dots, \beta_{5}). \end{aligned} \]
The transitions between states (30-day to 90-day and 90-day to 180-day) are also assigned Dirichlet distributions. In particular, for each state \(k\) at day 30 we place a Dirichlet prior distribution on the vector of probabilities for that state to the 90-day state \(j\), \(j = 0, \dots, 5\) under each treatment group separately: \[ \begin{aligned} &(p_{ctr,k,0}^{(30,90)}, \dots, p_{ctr,k,5}^{(30,90)} ) \sim \text{Dirichlet}(\alpha_{k,0}, \dots, \alpha_{k,5}), \quad k = 0, \dots, 5, \\ &(p_{trt,k,0}^{(30,90)}, \dots, p_{trt,k,5}^{(30,90)} ) \sim \text{Dirichlet}(\alpha_{k,0}, \dots, \alpha_{k,5}), \quad k = 0, \dots, 5. \end{aligned} \] Likewise, for each state \(j\) at 90 days we place a Dirichlet prior distribution on the vector of probabilities for that state to the 180-day state \(\ell\), \(\ell = 0, \dots, 5\) under each treatment group separately: \[ \begin{aligned} &(p_{ctr,j,0}^{(90, 180)}, \dots, p_{ctr,j,5}^{(90, 180)} ) \sim \text{Dirichlet}(\alpha_{j,0}, \dots, \alpha_{j,5}), \quad j = 0, \dots, 5, \\ &(p_{trt,j,0}^{(90, 180)}, \dots, p_{trt,j,5}^{(90, 180)} ) \sim \text{Dirichlet}(\alpha_{j,0}, \dots, \alpha_{j,5}), \quad j = 0, \dots, 5. \end{aligned} \]
These prior distributions are straightforward to update. The corresponding posteriors for the initial states are \[ \begin{aligned} &(p_{ctr,0}^{(30)}, \dots, p_{ctr,5}^{(30)} ) \mid (N_{ctr,0}^{(30)}, \dots, N_{ctr,5}^{(30)}) \sim \text{Dirichlet}(\beta_{0} + N_{ctr,0}^{(30)}, \dots, \beta_{5} + N_{ctr,5}^{(30)}), \\ &(p_{trt,0}^{(30)}, \dots, p_{trt,5}^{(30)} ) \mid (N_{trt,0}^{(30)}, \dots, N_{trt,5}^{(30)}) \sim \text{Dirichlet}(\beta_{0} + N_{trt,0}^{(30)}, \dots, \beta_{5} + N_{trt,5}^{(30)}), \end{aligned} \] where \(N_{trt,k}^{(30)}\) is the number of subjects in the treatment arm that had 30-day response value equal to \(k\). The posteriors for the 30 to 90 days transitions are \[ \begin{aligned} &(p_{ctr,k,0}^{(30,90)}, \dots, p_{ctr,k,5}^{(30,90)} ) \mid (N_{ctr,k,0}^{(30,90)}, \dots, N_{ctr,k,5}^{(30,90)}) \sim \text{Dirichlet}(\alpha_{k,0} + N_{ctr,k,0}^{(30,90)}, \dots, \alpha_{k,5} + N_{ctr,k,5}^{(30,90)}), \quad k = 0, \dots, 5, \\ &(p_{trt,k,0}^{(30,90)}, \dots, p_{trt,k,5}^{(30,90)} ) \mid (N_{trt,k,0}^{(30,90)}, \dots, N_{trt,k,5}^{(30,90)}) \sim \text{Dirichlet}(\alpha_{k,0} + N_{trt,k,0}^{(30,90)}, \dots, \alpha_{k,5} + N_{trt,k,5}^{(30,90)}), \quad k = 0, \dots, 5, \end{aligned} \] where \(N_{trt,k,j}^{(30,90)}\) is the number of subjects in the treatment arm that transitioned from 30-day state equal to \(k\) to 30-day state equal to \(j\). The posteriors for the 90 to 180 days transitions are \[ \begin{aligned} &(p_{ctr,j,0}^{(90, 180)}, \dots, p_{ctr,j,5}^{(90, 180)} ) \mid (N_{ctr,j,0}^{(90, 180)}, \dots, N_{ctr,j,5}^{(90, 180)}) \sim \text{Dirichlet}(\alpha_{j,0} + N_{ctr,j,0}^{(90, 180)}, \dots, \alpha_{j,5} + N_{ctr,j,5}^{(90, 180)}), \quad j = 0, \dots, 5, \\ &(p_{trt,j,0}^{(90, 180)}, \dots, p_{trt,j,5}^{(90, 180)} ) \mid (N_{trt,j,0}^{(90, 180)}, \dots, N_{trt,j,5}^{(90, 180)}) \sim \text{Dirichlet}(\alpha_{j,0} + N_{trt,j,0}^{(90, 180)}, \dots, \alpha_{j,5} + N_{trt,j,5}^{(90, 180)}), \quad j = 0, \dots, 5, \end{aligned} \] where \(N_{trt,j,\ell}^{(90,180)}\) is the number of subjects in the treatment arm that transitioned from 90-day state equal to \(j\) to 180-day state equal to \(\ell\).
3 Data Generation
To demonstrate how we can make inference on this model, we first need to simulate some data.
Let us load some useful functions and create the quantities corresponding to the underlying true model parameter.
Click to explore the functions in utils_OR.R.
installer <- function(package){
# install a package if it isn't installed
installed_packages <- package %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
install.packages(packages[!installed_packages])
}
}
# List of dependent packages
packages = c(
'MASS',
'ggalluvial',
'pammtools',
'dplyr',
'tidyr',
'stringr',
'purrr',
'ggplot2',
'latex2exp',
'RColorBrewer',
'rstan',
'bayesplot',
'coda',
'egg',
'parallel',
'forcats',
'MCMCpack',
'knitr',
'kableExtra'
)
# Install + source the packages
installer(packages)
invisible(lapply(packages, library, character.only = TRUE))
# Initialize global variables
if ('rstan' %in% packages) {
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores()) # use multiple cores if available, speeds up simulations
}
theme_set(theme_bw(base_size = 14))
options(dplyr.summarise.inform = FALSE)
# Auxiliary ---------------------------------------------------------------
#' Inverse logistic function
#' @param x: argument
#' @return inverse logistic applied to x
inv_logit = function(x) {
p = 1 / (1 + exp(-x))
p[p < 0] = exp(x) / (1 + exp(-1))
return(p)
}
#' Logistic function
#' @param p: argument
#' @return logistic applied to p
logit = function(p){
log(p / (1-p))
}
#' Probability density function of the logistic distribution
#' @param x: argument
#' @return pdf of the logistic applied to x
logistic_pdf = function(x) {
exp(-x) / (1 + exp(-x))^2
}
#' Rename the columns of a dataframe
#' @param df: data.frame whose column are to be renamed
#' @param names: names of the columns
#' @return data.frame with new column names
colnames_inplace = function(df, names){
colnames(df) = names
df
}
#' Calculate the class proportions in the treatment group for ordinal logistic regression
#' @param p_ctr: class proportions in the control group
#' @param OR: common odds ratio
#' @return class proportions in the treatment group
OR_transform = function(p_ctr, OR){
gamma = logit(head(cumsum(p_ctr), -1)) # transform to latent parameters (cutpoints)
theta = log(OR) # common Odds-ratio
p_trt = round(diff(c(0, inv_logit(theta + gamma), 1)), 2)
return (p_trt)
}
# Data generation ---------------------------------------------------------
#' Function to calculate the integrated accrual profile
#' @param t: time
#' @param peak_rate: peak accrual rate (weeks)
#' @param ramp_up: ramp up complete at week
#' @param t0: date of last enrolled patient (in weeks, since start of trial)
#' @return data.frame with new column names
Lambda = function(t, peak_rate, ramp_up, t0 = 0) {
if (ramp_up == 0){
peak_rate * (t + t0)
}
else {
0.5 * peak_rate/ramp_up * (t + t0)^2 * as.numeric((t + t0) < ramp_up) +
(peak_rate * ((t + t0) - ramp_up) + 0.5 * peak_rate * ramp_up) * as.numeric((t + t0) >= ramp_up)
}
}
#' Normalized integrated accrual profile
#' @param x: time
#' @param T_max: maximum follow up time
#' @param peak_rate: peak accrual rate (weeks)
#' @param ramp_up: ramp up complete at week
#' @param t0: date of last enrolled patient (in weeks, since start of trial)
#' @return data.frame with new column names
Ft = function(x, T_max, peak_rate, ramp_up, t0 = 0) {
(Lambda(x, peak_rate, ramp_up, t0) - Lambda(0, peak_rate, ramp_up, t0)) / (Lambda(T_max, peak_rate, ramp_up, t0) - Lambda(0, peak_rate, ramp_up, t0))
}
#' Inverse of the normalized integrated accrual profile
#' @param Ft: Normalized integrated accrual profile
#' @param u: uniform draw (to use the inverse cdf method)
#' @param T_max: maximum follow up time
#' @param peak_rate: peak accrual rate (weeks)
#' @param ramp_up: ramp up complete at week
#' @param t0: date of last enrolled patient (in weeks, since start of trial)
#' @return data.frame with new column names
Ftinv = function(Ft, u, T_max, peak_rate, ramp_up, t0 = 0) {
N = length(u)
a = rep(0, N)
b = rep(T_max, N)
binf = rep(0, N)
bsup = rep(0, N)
for (j in 1:20) {
idx_low = which(Ft((a + b) / 2, T_max, peak_rate, ramp_up, t0) <= u)
binf[idx_low] = (a[idx_low] + b[idx_low])/2
bsup[idx_low] = b[idx_low]
idx_upp = which(Ft((a + b) / 2, T_max, peak_rate, ramp_up, t0) >= u)
bsup[idx_upp] = (a[idx_upp] + b[idx_upp]) / 2
binf[idx_upp] = a[idx_upp]
a = binf
b = bsup
}
return ( (a + b) / 2 )
}
#' @param N: sample size
#' @param TRTPN: treatment arm indicator
#' @param p_ctr_true: probabilities under the control
#' @param p_trt_true: probabilities under the treatment
#' @param P_tr: transition probability matrix from Visit t to Visit (t-1)
#' @param T_max: follow up time after the last patient randomized
#' @param accr_profile: accrual profile
#' @param t0: date of last enrolled patient (in weeks, since start of trial)
#' @return
sample_long_ordinal_data = function(N, TRTPN, p_ctr_true, p_trt_true, P_tr, T_max, accr_profile, t0){
K = length(p_ctr_true)
# Sample times from start of accrual according to a Poisson Process
Ntimes = 5
n_max = rpois(1, Lambda(Ntimes*N/accr_profile$peak_rate, accr_profile$peak_rate, accr_profile$ramp_up, t0))
if (n_max < N){
stop('Cannot simulate n individuals with such small accrual!\n')
}
# Use the inverse cdf method to sample n_max accrual dates
dates_temp = Ftinv(Ft, runif(n_max), Ntimes*N/accr_profile$peak_rate, accr_profile$peak_rate, accr_profile$ramp_up, t0)
# Take only n accrual dates, sort them and shift them by the initial date
dates = t0 + sort(dates_temp)[1:N]
OTC30FL = RESP30 = OTC90FL = RESP90 = OTC180FL = RESP180 = array(NA, dim = N)
for (i in 1:N){
obs_window_i = dates[length(dates)] - dates[i] + T_max
if (obs_window_i < (180 + 21) / 7) { # if not past 180 days
OTC180FL[i] = 0
if (obs_window_i < (90 + 14) / 7) { # if not past 90 days
OTC90FL[i] = 0
if (obs_window_i < (30 + 14) / 7) { # if not past 30 days
OTC30FL[i] = 0
}
else { # between 30 and 90 days
OTC30FL[i] = 1
if (TRTPN[i] == 'Control') {
RESP30[i] = sample(0:(K - 1), 1, FALSE, p_ctr_true)
}
else if (TRTPN[i] == 'Treatment') {
RESP30[i] = sample(0:(K - 1), 1, FALSE, p_trt_true)
}
}
}
else { # between 90 and 180 days
OTC30FL[i] = OTC90FL[i] = 1
if (TRTPN[i] == 'Control') {
RESP90[i] = sample(0:(K - 1), 1, FALSE, p_ctr_true)
}
else if (TRTPN[i] == 'Treatment') {
RESP90[i] = sample(0:(K - 1), 1, FALSE, p_trt_true)
}
RESP30[i] = sample(0:(K - 1), 1, FALSE, P_tr[RESP90[i] + 1,])
}
}
else {
OTC30FL[i] = OTC90FL[i] = OTC180FL[i] = 1
if (TRTPN[i] == 'Control') {
RESP180[i] = sample(0:(K - 1), 1, FALSE, p_ctr_true)
}
else if (TRTPN[i] == 'Treatment') {
RESP180[i] = sample(0:(K - 1), 1, FALSE, p_trt_true)
}
RESP90[i] = sample(0:(K - 1), 1, FALSE, P_tr[RESP180[i] + 1,])
RESP30[i] = sample(0:(K - 1), 1, FALSE, P_tr[RESP90[i] + 1,])
}
}
data = cbind(dates, RESP30, RESP90, RESP180, OTC30FL, OTC90FL, OTC180FL) %>%
as_tibble %>%
mutate(SITEID = sample(paste0('0', 101:110), N, TRUE),
TRTPN = factor(key_trt[dose_label + 1]),
RANDDT = as.Date('2018-03-15') + 7 * dates,
CUTDT = max(RANDDT),
.before = dates) %>%
group_by(SITEID) %>%
mutate(SUBJID = paste0(SITEID, '-', str_pad(row_number(), 3, pad = "0")), .before = SITEID) %>%
ungroup() %>%
select(-dates)
return (data)
}
#' Function to plot the transition probability matrix
#' @param P_trans: list with three elements: p30 (vector), p90 (matrix), p180 (matrix)
#' @param digits: number of decimal digits to be used in the plot
#' @param title: title above the heatmaps
#' @param palette: color palette
#' @return plot object
plot_Ptrans = function(P_trans, digits = 1, title = NULL, palette = 'Oranges'){
K = length(P_trans$p30)
n_digits = if_else(digits == 0, "%s", paste0("%0.", digits, "f"))
# Initial state probabilities
p1 = P_trans$p30 %>%
as_tibble() %>%
mutate(rowname = as.character(0:(K - 1)),
colname = as.character(0), .before = 1) %>%
mutate(p_transitions = value / sum(value)) %>%
mutate(label = sprintf(n_digits, value)) %>%
ggplot(aes(colname, forcats::fct_rev(rowname), fill = p_transitions)) +
geom_tile() +
geom_text(aes(label = label)) +
coord_fixed() +
labs(x = '', y = 'Response at 30 days') +
scale_fill_distiller(palette = palette, direction = 1) +
theme(legend.position = 'none',
axis.text.x = element_blank(),
axis.ticks.x = element_blank()
)
# Transitions from 30 to 90
p2 = P_trans$p90 %>%
as_tibble() %>%
colnames_inplace(as.character(0:(K - 1))) %>%
mutate(rowname = as.character(0:(K - 1)), .before = 1) %>%
pivot_longer(
cols = -rowname,
names_to = "colname",
values_to = "transitions"
) %>%
group_by(rowname) %>%
mutate(p_transitions = transitions / sum(transitions),
p_transitions = if_else(is.na(p_transitions), 0, p_transitions)) %>%
ungroup() %>%
mutate(label = sprintf(n_digits, transitions)) %>%
ggplot(aes(colname, forcats::fct_rev(rowname), fill = p_transitions)) +
geom_tile() +
geom_text(aes(label = label)) +
coord_fixed() +
labs(x = 'Response at 90 days', y = 'Response at 30 days') +
scale_fill_distiller(palette = palette, direction = 1) +
theme(legend.position = 'none')
# Transitions from 90 to 180
p3 = P_trans$p180 %>%
as_tibble() %>%
colnames_inplace(as.character(0:(K - 1))) %>%
mutate(rowname = as.character(0:(K - 1)), .before = 1) %>%
pivot_longer(
cols = -rowname,
names_to = "colname",
values_to = "transitions"
) %>%
group_by(rowname) %>%
mutate(p_transitions = transitions / sum(transitions),
p_transitions = if_else(is.na(p_transitions), 0, p_transitions)) %>%
ungroup() %>%
mutate(label = sprintf(n_digits, transitions)) %>%
ggplot(aes(colname, forcats::fct_rev(rowname), fill = p_transitions)) +
geom_tile() +
geom_text(aes(label = label)) +
coord_fixed() +
labs(x = 'Response at 180 days', y = 'Response at 90 days') +
scale_fill_distiller(palette = palette, direction = 1) +
theme(legend.position = 'none')
plt = egg::ggarrange(p1, p2, p3,
ncol = 3,
widths = c(1, 6, 6),
top = title
)
return (plt)
}
# Model fit ---------------------------------------------------------------
#' Creates a dataframe with the posterior parameters for the Dirichlet
#' distribution used in the imputation model, for the two possible arms.
#' @param df: dataframe of study data
#' @param long_prior: prior for the Dirichlet distributions of the longitudinal model
#' @param key_trt: treatment arms labels
#' @return list of posterior hyperparameters for the Dirichlet distributions
get_imputation_pars = function(df, long_prior, key_trt){
# Number of outcome classes
K = length(long_prior$beta)
# Count how many transitions we have at 30, 90 and 180 days
PP30 = df %>%
filter(OTC30FL == 1) %>%
mutate(RESP30 = factor(RESP30, levels = 0:(K-1))) %>%
count(TRTPN, RESP30, .drop = F) %>%
filter(!is.na(RESP30))
PP90 = df %>%
filter(OTC90FL == 1) %>%
mutate(RESP30 = factor(RESP30, levels = 0:(K-1)),
RESP90 = factor(RESP90, levels = 0:(K-1))) %>%
count(TRTPN, RESP30, RESP90, .drop = F) %>%
filter(!is.na(RESP30), !is.na(RESP90))
PP180 = df %>%
filter(OTC180FL == 1) %>%
mutate(RESP90 = factor(RESP90, levels = 0:(K-1)),
RESP180 = factor(RESP180, levels = 0:(K-1))) %>%
count(TRTPN, RESP90, RESP180, .drop = F) %>%
filter(!is.na(RESP90), !is.na(RESP180))
# Calculate the updated parameters for the Dirichlet distributions
p30_ctr = PP30 %>% filter(TRTPN == key_trt[1]) %>% pull(n) + long_prior$beta
p30_trt = PP30 %>% filter(TRTPN == key_trt[2]) %>% pull(n) + long_prior$beta
p90_ctr = PP90 %>% filter(TRTPN == key_trt[1]) %>% pull(n) %>% matrix(K, K, byrow = T) + long_prior$alpha
p90_trt = PP90 %>% filter(TRTPN == key_trt[2]) %>% pull(n) %>% matrix(K, K, byrow = T) + long_prior$alpha
p180_ctr = PP180 %>% filter(TRTPN == key_trt[1]) %>% pull(n) %>% matrix(K, K, byrow = T) + long_prior$alpha
p180_trt = PP180 %>% filter(TRTPN == key_trt[2]) %>% pull(n) %>% matrix(K, K, byrow = T) + long_prior$alpha
return(list('p30_ctr' = p30_ctr,
'p30_trt' = p30_trt,
'p90_ctr' = p90_ctr,
'p90_trt' = p90_trt,
'p180_ctr' = p180_ctr,
'p180_trt' = p180_trt))
}
#' Main model function. At each iteration it imputes a dataset and performs the
#' relevant tests that declare success/futility for the clinical trial.
#' @param df: dataframe of study data
#' @param long_prior: prior for the Dirichlet distributions of the longitudinal model
#' @param n_max: maximum sample size
#' @param key_trt: treatment arms labels
#' @param iters: number of iterations
#' @param seed: random seed
#' @return list with elements:
#' 'PPn', a vector of flags for success/failure of the trials across simulations based on the rule PPn
#' 'PPmax', a vector of flags for success/failure of the trials across simulations based on the rule PPmax
#' 'OR_n', a vector with the observed ORs across simulations based on current sample size
#' 'OR_max', a vector with the observed ORs across simulations based on maximum sample size
#' 'pvals_n', a vector with the p-values across simulations based on maximum sample size
#' 'pvals_max', a vector with the p-values across simulations based on maximum sample size
fit_model = function(df, long_prior, n_max, key_trt, iters = 1000, seed = 12345){
# Prepare data structures
K = length(long_prior$beta)
PPn = array(NA, dim = iters)
PPmax = array(NA, dim = iters)
OR_n = array(NA, dim = iters)
OR_max = array(NA, dim = iters)
pvals_n = array(NA, dim = iters)
pvals_max = array(NA, dim = iters)
n_act = df %>%
pull(SUBJID) %>%
n_distinct()
ctr_act = df %>%
filter(TRTPN == key_trt[1]) %>%
nrow
trt_act = df %>%
filter(TRTPN == key_trt[2]) %>%
nrow
ctr_max = floor(n_max/3)
trt_max = n_max - ctr_max
n_ctr = ctr_max - ctr_act
n_trt = trt_max - trt_act
if ( (n_ctr <= 0) | (n_trt <= 0) ) {
stop("Number of patients in one of the arms is larger than maximum sample size\n")
}
# Get posterior parameters for the imputation model
imp_pars = df %>%
get_imputation_pars(long_prior, key_trt)
set.seed(seed)
# Sample all of the parameter vectors before the simulations loop
pi30_ctr = MCMCpack::rdirichlet(iters, imp_pars$p30_ctr)
pi30_trt = MCMCpack::rdirichlet(iters, imp_pars$p30_trt)
pi90_ctr = pi90_trt = pi180_ctr = pi180_trt = array(NA, dim = c(iters, K, K))
for (d in 1:K) {
pi90_ctr[,d,] = MCMCpack::rdirichlet(iters, imp_pars$p90_ctr[d,])
pi90_trt[,d,] = MCMCpack::rdirichlet(iters, imp_pars$p90_trt[d,])
pi180_ctr[,d,] = MCMCpack::rdirichlet(iters, imp_pars$p180_ctr[d,])
pi180_trt[,d,] = MCMCpack::rdirichlet(iters, imp_pars$p180_trt[d,])
}
# Calculate the relevant quantities for the current data under treatment and control
df_ctr = df %>%
filter(TRTPN == key_trt[1]) %>%
select(SUBJID, TRTPN, RESP30, RESP90, RESP180)
df_trt = df %>%
filter(TRTPN == key_trt[2]) %>%
select(SUBJID, TRTPN, RESP30, RESP90, RESP180)
# Some things to time the implementation
pb = txtProgressBar(min = 1, max = iters)
start = proc.time()
set.seed(seed)
my_seeds = sample(1:100000, iters, FALSE)
for (i in 1:iters){ # loop over simulations
set.seed(my_seeds[i])
### (1) Control group missing data imputation ###
pi30_temp = pi30_ctr[i,]
pi90_temp = pi90_ctr[i,,]
pi180_temp = pi180_ctr[i,,]
# Imputes the current dataset for the treatment group
df_imp_ctr = df_ctr %>%
add_row(
tibble(
SUBJID = paste0('imp_ctr_', 1:n_ctr),
TRTPN = key_trt[1]
)
) %>%
group_by(SUBJID) %>%
mutate(RESP30 = if_else(is.na(RESP30), as.numeric(sample(0:(K-1), 1, TRUE, pi30_temp)), RESP30),
RESP90 = if_else(is.na(RESP90), as.numeric(sample(0:(K-1), 1, TRUE, pi90_temp[RESP30 + 1,])), RESP90),
RESP180 = if_else(is.na(RESP180), as.numeric(sample(0:(K-1), 1, TRUE, pi180_temp[RESP90 + 1,])), RESP180)) %>%
ungroup()
### (2) Treatment group missing data imputation ###
pi30_temp = pi30_trt[i,]
pi90_temp = pi90_trt[i,,]
pi180_temp = pi180_trt[i,,]
# Imputes the current dataset for the treatment group
df_imp_trt = df_trt %>%
add_row(
tibble(
SUBJID = paste0('imp_trt_', 1:n_trt),
TRTPN = key_trt[2]
)
) %>%
group_by(SUBJID) %>%
mutate(RESP30 = if_else(is.na(RESP30), as.numeric(sample(0:(K-1), 1, TRUE, pi30_temp)), RESP30),
RESP90 = if_else(is.na(RESP90), as.numeric(sample(0:(K-1), 1, TRUE, pi90_temp[RESP30 + 1,])), RESP90),
RESP180 = if_else(is.na(RESP180), as.numeric(sample(0:(K-1), 1, TRUE, pi180_temp[RESP90 + 1,])), RESP180)) %>%
ungroup()
df_imp = rbind(df_imp_trt, df_imp_ctr)
### (3) Perform the ordinal regression test ###
# Current sample size
data_fit = df_imp %>%
filter(!str_detect(SUBJID, 'imp_')) %>%
mutate(RESP180 = factor(RESP180))
fit_freq = MASS::polr(RESP180 ~ TRTPN, data = data_fit, Hess = TRUE)
fit_summ = summary(fit_freq)
obs_OR = exp(- fit_summ$coefficients[1,"Value"])
OR_n[i] = obs_OR
t_val = - fit_summ$coefficients[1,"Value"] / fit_summ$coefficients[1,"Std. Error"]
dof = fit_freq$df.residual
pval_n = 1 - pt(q = t_val, df = dof)
pvals_n[i] = pval_n
# Maximum sample size
data_fit = df_imp %>%
mutate(RESP180 = factor(RESP180))
fit_freq = MASS::polr(RESP180 ~ TRTPN, data = data_fit, Hess = TRUE)
fit_summ = summary(fit_freq)
obs_OR = exp(- fit_summ$coefficients[1,"Value"])
OR_max[i] = obs_OR
t_val = - fit_summ$coefficients[1,"Value"] / fit_summ$coefficients[1,"Std. Error"]
dof = fit_freq$df.residual
pval_max = 1 - pt(q = t_val, df = dof)
pvals_max[i] = pval_max
# Calculate rules to declare success/futility at the final
PPn[i] = if_else(pval_n < 0.02, 1, 0)
PPmax[i] = if_else(pval_max < 0.02, 1, 0)
setTxtProgressBar(pb, i)
}
cat(sprintf('\n Imputation complete in %.1f minutes', (proc.time() - start)[3]/ 60))
return(list(
PPn = PPn,
PPmax = PPmax,
OR_n = OR_n,
OR_max = OR_max,
pvals_n = pvals_n,
pvals_max = pvals_max
))
}We define the true control ordinal probabilities p_ctr_true and the true odds ratio OR_true.
These two quantities imply the true treatment ordinal probabilities p_trt_true that can be calculated with the function OR_transform().
source('./utils_OR.R')
cols = brewer.pal(8, 'Dark2')
key_trt = c('Control', 'Treatment')
data_seed = 4321
weeks_to_months = 6/26
p_ctr_true = c(0.07, 0.2, 0.28, 0.2, 0.15, 0.1) # true control ordinal probabilities
K = length(p_ctr_true) # number of ordinal categories (from 0 to K-1)
OR_true = 1.8 # true odds ratio
p_trt_true = OR_transform(p_ctr_true, OR_true) # true implied treatment ordinal probabilitiesWe first have to generate the allocation to each arm. In our example, we use an allocation 1:1 in blocks of 2 so that at the analysis sample size we have a perfectly balanced allocation.
# Generate the data
set.seed(data_seed)
N = 300
dose_label = sample(c(rep(0, N/2), rep(1, N/2)))
TRTPN = factor(key_trt[dose_label + 1], levels = key_trt)We can also setup the accrual parameters. For example, consider a peak rate accrual of 12 patients per month with a 24-month ramp-up period.
# Accrual profile
accr_profile = NULL
accr_profile$peak_rate = 12 * weeks_to_months # in weeks
accr_profile$ramp_up = 24 / weeks_to_months # in weeksTo generate the data, we will first generate the 180-day outcomes using the class probabilities that we just defined.
Afterwards, we generate previous visits backwards according to the transition probabilities defined below in P_tr.
This matrix has on each row the probabilities of the response class at the previous visits given the response at the current visit.
# Transition probabilities
P_tr = matrix(c(0.7, 0.3, 0, 0, 0, 0,
0, 0.2, 0.4, 0.4, 0, 0,
0, 0, 0.34, 0.66, 0, 0,
0, 0, 0.22, 0.67, 0.11, 0,
0, 0, 0, 0.71, 0.29, 0,
0, 0, 0, 0, 0.3, 0.7),
K, K, byrow = T)
# Show P_tr in output
rownames(P_tr) = as.character(0:(K-1))
P_tr %>%
knitr::kable(booktabs = TRUE,
col.names = as.character(0:(K-1)),
caption = 'Probability of response values at previous visit
given score at current visit.') %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"))| 0 | 1 | 2 | 3 | 4 | 5 | |
|---|---|---|---|---|---|---|
| 0 | 0.7 | 0.3 | 0.00 | 0.00 | 0.00 | 0.0 |
| 1 | 0.0 | 0.2 | 0.40 | 0.40 | 0.00 | 0.0 |
| 2 | 0.0 | 0.0 | 0.34 | 0.66 | 0.00 | 0.0 |
| 3 | 0.0 | 0.0 | 0.22 | 0.67 | 0.11 | 0.0 |
| 4 | 0.0 | 0.0 | 0.00 | 0.71 | 0.29 | 0.0 |
| 5 | 0.0 | 0.0 | 0.00 | 0.00 | 0.30 | 0.7 |
Let us generate the data with the function sample_long_ordinal_data().
set.seed(data_seed)
data = sample_long_ordinal_data(N = N,
TRTPN,
p_ctr_true = p_ctr_true,
p_trt_true = p_trt_true,
P_tr,
T_max = 0,
accr_profile = accr_profile,
t0 = 0)We can visualize the data we just generated with a series of plots. For example, a stacked barplot of the response categories under both treatment arms can be useful to see the empirical cumulative distribution of the response at the final visit.
data %>%
filter(OTC180FL == 1) %>%
group_by(TRTPN, RESP180) %>%
summarise(n = n()) %>%
mutate(freq = n / sum(n)) %>%
right_join(list(
TRTPN = levels(TRTPN),
RESP180 = 0:(K-1)) %>%
cross_df(),
by = c('TRTPN', 'RESP180')) %>%
mutate(TRTPN = factor(TRTPN),
n = ifelse(is.na(n), 0, n),
freq = ifelse(is.na(freq), 0, freq)) %>%
arrange(TRTPN, RESP180) %>%
group_by(TRTPN) %>%
mutate(freq_sum = cumsum(freq)) %>%
ggplot() +
geom_bar(aes(x = freq, y = TRTPN, fill = factor(RESP180, levels = 0:(K-1))),
position = position_stack(reverse = TRUE), stat = 'identity') +
labs(x = 'Proportion of subjects',
y = '',
title = 'Frequency of the response at 180 days') +
scale_fill_brewer(name = 'Response', palette = "RdBu", direction=-1) +
theme(legend.position = "top") +
guides(fill = guide_legend(nrow = 1))
Equivalently, a simple barplot shows the empirical distribution of the response classes at the final visit.
data %>%
filter(OTC180FL == 1) %>%
group_by(TRTPN, RESP180) %>%
summarise(n = n()) %>%
mutate(freq = n / sum(n)) %>%
right_join(list(
TRTPN = levels(TRTPN),
RESP180 = 0:(K-1)) %>%
cross_df(),
by = c('TRTPN','RESP180')) %>%
mutate(TRTPN = factor(TRTPN),
n = ifelse(is.na(n), 0, n),
freq = ifelse(is.na(freq), 0, freq)) %>%
arrange(TRTPN, RESP180) %>%
group_by(TRTPN) %>%
mutate(freq_sum = cumsum(freq),
imp_month_plot = ifelse(TRTPN == key_trt[1], RESP180, RESP180)) %>%
ggplot() +
geom_step(aes(x = imp_month_plot, y = freq_sum, col = TRTPN), size = 1) +
geom_bar(aes(fill = TRTPN, col = TRTPN, y = freq, x = RESP180),
width = 0.6, position = position_dodge(width = 0.7), stat = "identity") +
labs(x = 'Proportion of subjects',
y = '',
title = 'Frequency of the response at 180 days') +
scale_x_continuous(breaks = 0:(K-1)) +
scale_colour_brewer(name = '', palette = 'Set2') +
scale_fill_brewer(name = '', palette = 'Set2') +
theme(legend.position = "top")
We can also count the number of transitions from each pairs of response variables and show them for each group separately. This shows what transitions are most common from each of the observed states. The two figures below are color shaded according to the probability of transitioning from the current state to the possible future states. Note that the color shading is normalized in each row, i.e., it denotes the probability of transitioning to states given each fixed current state.
# Set prior counts to 0 so that only the empirical transitions are
# counted
pr_counts = list(beta = rep(0, K),
alpha = matrix(c(0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0),
K, K, byrow = T))
# Updated parameters for the longitudinal model
imp_pars = data %>%
select(TRTPN,
OTC30FL, RESP30,
OTC90FL, RESP90,
OTC180FL, RESP180) %>%
get_imputation_pars(pr_counts, key_trt)
# Plot the number of transitions under control group
plt_ctr = plot_Ptrans(list(p30 = imp_pars$p30_ctr,
p90 = imp_pars$p90_ctr,
p180 = imp_pars$p180_ctr),
digits = 0,
palette = 'Greens')
# Plot the number of transitions under treatment group
plt_trt = plot_Ptrans(list(p30 = imp_pars$p30_trt,
p90 = imp_pars$p90_trt,
p180 = imp_pars$p180_trt),
digits = 0)
A more comprehensive visualization for ordinal outcomes that shows simultaneously the distribution of the responses at each visit and the transitions between visits is the alluvial plot below.
data %>%
group_by(TRTPN) %>%
count(RESP30, RESP90, RESP180) %>%
mutate(id = row_number()) %>%
pivot_longer(cols = contains('RESP')) %>%
mutate(value = factor(value, levels = c(NA, as.character(0:(K-1))), exclude = NULL),
name = factor(name, levels = c('RESP30', 'RESP90', 'RESP180'))) %>%
group_by(name) %>%
ggplot(aes(x = name, y = n, stratum = value, fill = value, alluvium = id)) +
ggalluvial::geom_stratum(alpha = 1) +
ggalluvial::geom_flow() +
facet_wrap(~ TRTPN) +
scale_fill_brewer(name = 'Response', breaks = c(NA, as.character(0:(K-1))), labels = c('Missing', as.character(0:(K-1))), palette = "RdBu", direction = -1, na.value = 'grey75') +
scale_x_discrete(name = '', labels = c('30 Days', '90 Days', '180 Days'), expand = c(0.1, 0.1)) +
scale_y_continuous(name = 'Proportion of patients', expand = c(0.01, 0.01))
We can also inspect the accrual of the trial that we just simulated.
data %>%
mutate(TIMESINCERAND = as.numeric(difftime(RANDDT, min(RANDDT), units = 'weeks'))) %>%
arrange(TIMESINCERAND) %>%
mutate(Subject = row_number(),
EXPACCR = Lambda(t = TIMESINCERAND, peak_rate = accr_profile$peak_rate,
ramp_up = accr_profile$ramp_up)) %>%
select(TIMESINCERAND, Subject, EXPACCR) %>%
ggplot() +
geom_step(aes(x = TIMESINCERAND, y = Subject, col = 'observed'), size = 1) +
geom_line(aes(x = TIMESINCERAND, y = EXPACCR,
col = sprintf('%.2f per week', accr_profile$peak_rate)),
size = 1) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
labs(x = 'Weeks since start of accrual',
y = 'Cumulative number of subjects',
title = 'Observed vs simulated accrual') +
scale_color_brewer(name = '', palette = 'Set1') +
theme(legend.position = "top")
4 Model Fitting
4.1 Frequentist Model Fitting
The first modeling strategy to fit the data that we just simulated is a frequentist ordinal logistic regression on the 180-day response values. This model uses information on completers only and excludes patients that are still in follow-up.
The ordinal regression can be fit using the MASS package in R.
# Fit the frequentist version
fit_freq = MASS::polr(factor(RESP180) ~ TRTPN, data = data, Hess = TRUE)
fit_summ = summary(fit_freq)Be aware that the polr() function uses
\[\text{logit} \{P(Y \leq k \mid x)\} = c_{k} - x\theta,\]
which is a slightly different parametrization from the one described in Section 2.1.
Thus, to recover the odds ratio we need to flip the sign of the coefficient returned to the polr() function, as done below.
# Compare model parameters with true values
log(OR_true); - fit_summ$coefficients[1,1][1] 0.5877867[1] 0.4718522logit(head(cumsum(p_ctr_true), -1)); as.numeric(fit_summ$coefficients[-1,1])[1] -2.5866893 -0.9946226 0.2006707 1.0986123 2.1972246[1] -2.2464667 -0.8908499 0.2816976 1.2053048 2.2016640As one can see, the model seems to be close to the true values for both the odds ratio and for the latent cutpoints.
Suppose that for a trial success is demonstrated if the one-sided p-value corresponding to the test that the common odds ratio (OR), \(e^{\theta}\), is greater than \(1\) is less than or equal to \(0.02\). We can calculate such p-value with the code below.
obs_OR = exp(- fit_summ$coefficients[1,1]) # observed OR
t_val = - fit_summ$coefficients[1,1] / fit_summ$coefficients[1,2] # t value for the test
dof = fit_freq$df.residual # degrees of freedom
pval = 1 - pt(q = t_val, df = dof)
pval; pval < 0.02[1] 0.02429251[1] FALSE4.2 Bayesian Model Fitting
An alternative model is a Bayesian ordinal logistic regression on the 180-day response values. This model also uses information on completers only and excludes patients that are still in follow-up.
4.2.1 Prior Elicitation
We need to set priors on all parameters of interest, i.e., the internal cutpoints and the odds ratio.
Unless response classes are sparse, i.e., there are empty classes, a uniform prior on the internal cutpoints suffices.
We can also use a uniform prior on the log odds ratio.
Priors that can be used for the longitudinal model are discussed in Section 7.
4.2.2 Posterior Sampling
The STAN model is the following:
writeLines(readLines("./ord_logistic.stan"))data {
int<lower=1> N; // Number of observations
int<lower=1> D; // Max of ordinal categories
int<lower=0, upper=1> X[N];
int<lower=1, upper=D> y[N]; // Observed ordinals {0, ..., K}
}
parameters {
real theta; // Latent effect
ordered[D - 1] c; // (Internal) cut points
}
model {
// Prior model
theta ~ normal(0, 2); // prior on the OR parameter
// Implicit uniform priors on the cutpoints
// Observational model
for (i in 1:N) {
y[i] ~ ordered_logistic(0.0, c + X[i] * theta);
}
}The function ordered_logistic_lpmf(0, c + X[i] * theta) is equivalent to the class probabilities \(p_{k} = \Pi(- c_{k−1} - \theta) - \Pi(- c_{k} - \theta)\).
We can now run the model and put the results (MCMC chains) into a matrix whose columns correspond to the different model parameters \(\theta, c_{1}, \dots, c_{K-1}\).
# Choose the seed
model_seed = 12345
# STAN data
stan_data <- list(N = nrow(data %>% filter(OTC180FL == 1)),
D = K,
X = as.numeric(data %>% filter(OTC180FL == 1) %>% pull(TRTPN)) - 1,
y = data %>% filter(OTC180FL == 1) %>% pull(RESP180) + 1
)
# Run STAN model
n_iters = 5000
burnin = 2500
n_thin = 2
n_chains = 4
samp_size = n_chains * (n_iters - burnin) / n_thin
set.seed(model_seed)
fit <- stan(file = 'ord_logistic.stan',
data = stan_data,
chains = n_chains,
warmup = burnin,
iter = n_iters,
thin = n_thin,
seed = model_seed) %>%
extract()
# Extract model output
samples = fit$theta %>%
cbind(fit$c) %>%
colnames_inplace(c('theta', sprintf('c[%s]', 1:(K - 1))))4.2.3 MCMC Diagnostics
The traceplots for the model parameters do not show any concerning behavior.
# Traceplots
mcmc_trace(samples, facet_args = list(labeller = ggplot2::label_parsed,
ncol = 2))
The total sample size of the chains was 5000, and the effective sample sizes are very close to their optimal values.
# Effective sample size
print(round(coda::effectiveSize(samples)))theta c[1] c[2] c[3] c[4] c[5]
4496 4750 4793 5000 5000 4779 5 Posterior Checks
Once we have the MCMC draws, we can calculate any functional of interest, e.g. posterior mean, median, quantiles, etc., and check that the estimated parameters are similar to the ground truth.
log(OR_true); mean(samples[,1])[1] 0.5877867[1] 0.4673618logit(head(cumsum(p_ctr_true), -1)); as.numeric(colMeans(samples[,2:K]))[1] -2.5866893 -0.9946226 0.2006707 1.0986123 2.1972246[1] -2.2793876 -0.8970090 0.2893855 1.2279495 2.2607256We can also illustrate the posterior distributions of the cutpoints.
# Show the posterior distributions of the cutpoints
samples %>%
as_tibble %>%
select(contains('c[')) %>%
pivot_longer(cols = everything(), names_to = 'cutpoint') %>%
ggplot() +
geom_histogram(aes(x = value, y = ..density.., fill = cutpoint),
col = 'white', alpha = 0.4) +
scale_x_continuous(name = TeX('$c_{j}$')) +
scale_fill_brewer(name = '', palette = 'Dark2', labels = sprintf('j = %i', 1:(K-1))) +
theme(legend.position = 'top') +
guides(fill = guide_legend(nrow = 1))
An effective way to show the predictive inference from the model fit is to show the posterior distributions of the class probabilities under control and treatment group.
# Show the posterior distribution of the class probabilities under the control group
p_ctr_post = cbind(rep(0, samp_size), inv_logit(samples[,2:K]), rep(1, samp_size)) %>%
apply(1, diff) %>%
t() %>%
colnames_inplace(sprintf('p_%s', 0:(K - 1)))
p_ctr_est = p_ctr_post %>%
as_tibble %>%
pivot_longer(cols = everything(), names_to = 'param') %>%
group_by(param) %>%
summarize(mean = mean(value),
low = quantile(value, 0.025),
upp = quantile(value, 0.975)) %>%
mutate(true = p_ctr_true,
param = as.numeric(str_remove(param, 'p_')))
p_ctr_est %>%
rbind(tail(p_ctr_est, 1) %>% mutate(param = param + 1)) %>%
ggplot() +
geom_step(aes(x = param - 0.5, y = mean), col = brewer.pal(8, 'Dark2')[1], size = 1) +
pammtools::geom_stepribbon(aes(x = param - 0.5, ymin = low, ymax = upp), fill = brewer.pal(8, 'Dark2')[1], alpha = 0.3) +
geom_step(aes(x = param - 0.5, y = true)) +
scale_y_continuous(name = '', expand = c(0, 0), limits = c(0, max(p_ctr_est$upp) + 0.1)) +
scale_x_continuous(name = '', breaks = 0:(K-1), labels = TeX(sprintf("$p_{\\%s}$", 0:(K-1))), expand = c(0, 0))
# Show the posterior distribution of the class probabilities under the treatment group
p_trt_post = cbind(rep(0, samp_size), inv_logit(samples[,2:K] + samples[,1]), rep(1, samp_size)) %>%
apply(1, diff) %>%
t() %>%
colnames_inplace(sprintf('p_%s', 0:(K - 1)))
p_trt_est = p_trt_post %>%
as_tibble %>%
pivot_longer(cols = everything(), names_to = 'param') %>%
group_by(param) %>%
summarize(mean = mean(value),
low = quantile(value, 0.025),
upp = quantile(value, 0.975)) %>%
mutate(true = p_trt_true,
param = as.numeric(str_remove(param, 'p_')))
p_trt_est %>%
rbind(tail(p_trt_est, 1) %>% mutate(param = param + 1)) %>%
ggplot() +
geom_step(aes(x = param - 0.5, y = mean), col = brewer.pal(8, 'Dark2')[2], size = 1) +
pammtools::geom_stepribbon(aes(x = param - 0.5, ymin = low, ymax = upp), fill = brewer.pal(8, 'Dark2')[2], alpha = 0.3) +
geom_step(aes(x = param - 0.5, y = true)) +
scale_y_continuous(name = '', expand = c(0, 0), limits = c(0, max(p_trt_est$upp) + 0.1)) +
scale_x_continuous(name = '', breaks = 0:(K-1), labels = TeX(sprintf("$p_{\\%s}$", 0:(K-1))), expand = c(0, 0))
The black lines represent the underlying truth, the colored lines and shaded areas are the mean and \(95\%\) credible interval from the Bayesian model.
6 Model Diagnostics
We propose here two graphic diagnostics to check the proportional odds assumption of the model. The first one involves fitting repeated logistic regressions to all possible dichotomizations of the ordinal categories and verify that the odds ratios are similar to each other.
# Fit the frequentist PO model
data_fit = data %>%
mutate(RESP180 = factor(RESP180))
fit_freq = MASS::polr(RESP180 ~ TRTPN, data = data_fit, Hess = TRUE)
fit_summ = summary(fit_freq)
obs_OR = exp(- fit_summ$coefficients[1,"Value"])
LI_OR = exp(- fit_summ$coefficients[1,"Value"] - 1.96 * fit_summ$coefficients[1,"Std. Error"])
UI_OR = exp(- fit_summ$coefficients[1,"Value"] + 1.96 * fit_summ$coefficients[1,"Std. Error"])
binned_OR = tibble('method' = 'Overall PO',
'mean' = obs_OR,
'li' = LI_OR,
'ui' = UI_OR)
# Fit the dichotomized logistic regressions
for (i in 0:(K-2)){
data_binn = data %>%
mutate(RESP180 = factor(if_else(RESP180 <= i, 0, 1)))
fit_binn = glm(RESP180 ~ TRTPN, data = data_binn, family = "binomial")
fit_summ_binn = summary(fit_binn)
obs_OR = exp(- fit_summ_binn$coefficients[2,"Estimate"])
LI_OR = exp(- fit_summ_binn$coefficients[2,"Estimate"] - 1.96 * fit_summ_binn$coefficients[2,"Std. Error"])
UI_OR = exp(- fit_summ_binn$coefficients[2,"Estimate"] + 1.96 * fit_summ_binn$coefficients[2,"Std. Error"])
binned_OR = binned_OR %>%
add_row('method' = paste0('Response >= ', i + 1, ' vs others'),
'mean' = obs_OR,
'li' = LI_OR,
'ui' = UI_OR)
}
# Overall odds ratio estimate from the PO model and individual
# dichotomized odds ratio estimates
binned_OR %>%
mutate(method = factor(method, levels = c(sprintf('Response >= %s vs others', (K-1):1), 'Overall PO'))) %>%
ggplot() +
geom_rect(aes(ymin = as.numeric(binned_OR[1,'li']), ymax = as.numeric(binned_OR[1,'ui']), xmin = -Inf, xmax = +Inf), alpha = 0.05) +
geom_linerange(aes(ymin = li, ymax = ui, x = method), lwd = 1) +
geom_point(aes(x = method, y = mean, col = "Estimated OR"), size = 3) +
geom_hline(yintercept = 1, lty = 5) +
geom_hline(yintercept = as.numeric(binned_OR[1,'mean']), lty = 3) +
coord_flip() +
scale_color_manual(name = '', values = brewer.pal(8, 'Dark2')[2]) +
scale_x_discrete(name = '') +
scale_y_continuous(name = "Odds Ratio (log scale)", trans='log2') +
theme(legend.position = 'none')
We can see here that the odds ratio are all similar to each other (e.g., within each other’s confidence intervals). Using the binary logistic models, moreover, we can see that small effects of the treatment are observed for each possible dichotomization of the ordinal outcome. However, these effects are usually not significant. Performing ordinal logistic regression, we can produce a common odds ratio, which has a narrower confidence interval. Thus, in this case this method has greater power to detect a significant effect.
We can also inspect the log-odds (logit transformed cumulative probabilities) and verify that they differ by a constant across the two treatment groups. The proportional odds model (dashed line), in fact, assumes that there is a constant vertical shift between the log-odds.
# Log-odds (logit transformed cumulative probabilities).
# The model assumes a constant vertical shift between the two groups
data_fit %>%
group_by(TRTPN, RESP180) %>%
summarize(n = n()) %>%
complete(RESP180) %>%
mutate(n = if_else(is.na(n), as.integer(0), n)) %>%
filter(!is.na(RESP180)) %>%
mutate(emp_p = n / sum(n),
cumul_emp_p = cumsum(emp_p)) %>%
ungroup() %>%
mutate(cumul_fitted_p = c(as.numeric(inv_logit(fit_summ$coefficients[-1,'Value'])),
1,
inv_logit(- fit_summ$coefficients[1,'Value'] + as.numeric(fit_summ$coefficients[-1,'Value'])),
1)) %>%
ggplot() +
geom_step(aes(x = as.numeric(RESP180), y = logit(cumul_emp_p), col = TRTPN), size = 1) +
geom_step(aes(x = as.numeric(RESP180), y = logit(cumul_fitted_p), col = TRTPN), lty = 2, size = 1) +
scale_x_continuous(name = 'Response at 180 days', breaks = 0:(K-1), labels = levels(data_fit$RESP180)) +
scale_y_continuous(name = 'logit(cumulative probabilities) \n [log-odds]') +
scale_color_brewer(name = '', palette = 'Dark2') +
theme(legend.position = 'top')
As we can see in this graph (and as we knew, since we generated the data), the proportional odds assumption is verified in the empirical data (solid lines).
7 Predictive Probabilities
Let’s now use a model that includes information from the previous visits of non-completers.
As before, suppose that we designed a trial with maximum sample size \(n_{max} = 500\) where final success is demonstrated if the one-sided p-value corresponding to the hypothesis that the common odds ratio (OR), \(e^{\theta}\), is greater than \(1\) is less than or equal to \(0.02\). This threshold value is chosen based on the early interim analyses to stop for predicted success to maintain an overall one-sided type I error rate of \(0.025\).
At each interim, we only observe partial data, i.e., not all subjects are completers.
For instance, in the data we simulated we have \(n = 300\) randomized subjects and not all of them are completers.
We can use predictive probabilities of success to determine if the trial should stop for early success or futility at the pre-planned interim analyses.
The predictive probability of success at the current sample size (\(PP_{n}\)) combines the observed completers with the uncertainty in the currently enrolled patients that have not completed yet.
The predictive probability of success at the maximum sample size (\(PP_{max}\)) combines the observed completers with the uncertainty in the currently enrolled patients that have not completed yet and the uncertainty in the future enrolled patients.
The predictive probabilities are calculated based on the longitudinal model described in Section 2.2 by the following procedure:
- For subjects without 30 day data, whether they have been enrolled (current sample size) or not (future patients): simulate 30-day response from posterior distribution of the response at 30 days
- For subjects without 90 day data: simulate 90-day response from posterior distribution of the response at 90 days given the observed or above predicted 30-day data and treatment arm
- For subjects without 180 day data: simulate 180-day response from posterior distribution of the response at 180 days given the observed or above predicted 90-day data and treatment arm
- Determine if the if the completed (observed and predicted) data at the current sample size and if the completed (observed and predicted) data at the maximum sample size lead to a successful trial (p-value corresponding to the test that \(\text{OR} > 1\) is less than or equal to \(0.02\))
- Repeat steps 1-4 \(1000\) times and report the proportion of times that the completed data at current and maximum sample size result in a successful trial, i.e., \(PP_{n}\) and \(PP_{max}\), respectively
We will use the function fit_model() illustrated below to calculate the predictive probability of success at the current sample size and at the maximum sample size.
print(fit_model)function (df, long_prior, n_max, key_trt, iters = 1000, seed = 12345)
{
K = length(long_prior$beta)
PPn = array(NA, dim = iters)
PPmax = array(NA, dim = iters)
OR_n = array(NA, dim = iters)
OR_max = array(NA, dim = iters)
pvals_n = array(NA, dim = iters)
pvals_max = array(NA, dim = iters)
n_act = df %>% pull(SUBJID) %>% n_distinct()
ctr_act = df %>% filter(TRTPN == key_trt[1]) %>% nrow
trt_act = df %>% filter(TRTPN == key_trt[2]) %>% nrow
ctr_max = floor(n_max/3)
trt_max = n_max - ctr_max
n_ctr = ctr_max - ctr_act
n_trt = trt_max - trt_act
if ((n_ctr <= 0) | (n_trt <= 0)) {
stop("Number of patients in one of the arms is larger than maximum sample size\n")
}
imp_pars = df %>% get_imputation_pars(long_prior, key_trt)
set.seed(seed)
pi30_ctr = MCMCpack::rdirichlet(iters, imp_pars$p30_ctr)
pi30_trt = MCMCpack::rdirichlet(iters, imp_pars$p30_trt)
pi90_ctr = pi90_trt = pi180_ctr = pi180_trt = array(NA, dim = c(iters,
K, K))
for (d in 1:K) {
pi90_ctr[, d, ] = MCMCpack::rdirichlet(iters, imp_pars$p90_ctr[d,
])
pi90_trt[, d, ] = MCMCpack::rdirichlet(iters, imp_pars$p90_trt[d,
])
pi180_ctr[, d, ] = MCMCpack::rdirichlet(iters, imp_pars$p180_ctr[d,
])
pi180_trt[, d, ] = MCMCpack::rdirichlet(iters, imp_pars$p180_trt[d,
])
}
df_ctr = df %>% filter(TRTPN == key_trt[1]) %>% select(SUBJID,
TRTPN, RESP30, RESP90, RESP180)
df_trt = df %>% filter(TRTPN == key_trt[2]) %>% select(SUBJID,
TRTPN, RESP30, RESP90, RESP180)
pb = txtProgressBar(min = 1, max = iters)
start = proc.time()
set.seed(seed)
my_seeds = sample(1:100000, iters, FALSE)
for (i in 1:iters) {
set.seed(my_seeds[i])
pi30_temp = pi30_ctr[i, ]
pi90_temp = pi90_ctr[i, , ]
pi180_temp = pi180_ctr[i, , ]
df_imp_ctr = df_ctr %>% add_row(tibble(SUBJID = paste0("imp_ctr_",
1:n_ctr), TRTPN = key_trt[1])) %>% group_by(SUBJID) %>%
mutate(RESP30 = if_else(is.na(RESP30), as.numeric(sample(0:(K -
1), 1, TRUE, pi30_temp)), RESP30), RESP90 = if_else(is.na(RESP90),
as.numeric(sample(0:(K - 1), 1, TRUE, pi90_temp[RESP30 +
1, ])), RESP90), RESP180 = if_else(is.na(RESP180),
as.numeric(sample(0:(K - 1), 1, TRUE, pi180_temp[RESP90 +
1, ])), RESP180)) %>% ungroup()
pi30_temp = pi30_trt[i, ]
pi90_temp = pi90_trt[i, , ]
pi180_temp = pi180_trt[i, , ]
df_imp_trt = df_trt %>% add_row(tibble(SUBJID = paste0("imp_trt_",
1:n_trt), TRTPN = key_trt[2])) %>% group_by(SUBJID) %>%
mutate(RESP30 = if_else(is.na(RESP30), as.numeric(sample(0:(K -
1), 1, TRUE, pi30_temp)), RESP30), RESP90 = if_else(is.na(RESP90),
as.numeric(sample(0:(K - 1), 1, TRUE, pi90_temp[RESP30 +
1, ])), RESP90), RESP180 = if_else(is.na(RESP180),
as.numeric(sample(0:(K - 1), 1, TRUE, pi180_temp[RESP90 +
1, ])), RESP180)) %>% ungroup()
df_imp = rbind(df_imp_trt, df_imp_ctr)
data_fit = df_imp %>% filter(!str_detect(SUBJID, "imp_")) %>%
mutate(RESP180 = factor(RESP180))
fit_freq = MASS::polr(RESP180 ~ TRTPN, data = data_fit,
Hess = TRUE)
fit_summ = summary(fit_freq)
obs_OR = exp(-fit_summ$coefficients[1, "Value"])
OR_n[i] = obs_OR
t_val = -fit_summ$coefficients[1, "Value"]/fit_summ$coefficients[1,
"Std. Error"]
dof = fit_freq$df.residual
pval_n = 1 - pt(q = t_val, df = dof)
pvals_n[i] = pval_n
data_fit = df_imp %>% mutate(RESP180 = factor(RESP180))
fit_freq = MASS::polr(RESP180 ~ TRTPN, data = data_fit,
Hess = TRUE)
fit_summ = summary(fit_freq)
obs_OR = exp(-fit_summ$coefficients[1, "Value"])
OR_max[i] = obs_OR
t_val = -fit_summ$coefficients[1, "Value"]/fit_summ$coefficients[1,
"Std. Error"]
dof = fit_freq$df.residual
pval_max = 1 - pt(q = t_val, df = dof)
pvals_max[i] = pval_max
PPn[i] = if_else(pval_n < 0.02, 1, 0)
PPmax[i] = if_else(pval_max < 0.02, 1, 0)
setTxtProgressBar(pb, i)
}
cat(sprintf("\n Imputation complete in %.1f minutes", (proc.time() -
start)[3]/60))
return(list(PPn = PPn, PPmax = PPmax, OR_n = OR_n, OR_max = OR_max,
pvals_n = pvals_n, pvals_max = pvals_max))
}The only remaining thing to do is to discuss the prior hyperparameters for the longitudinal model used in the imputation for the predictive probability calculations.
For the prior on the 30-day responses, we choose \(\beta_{k} = 1\) for \(k = 0, \dots, 5\) and we use the same prior by treatment arm.
Note, however, that the probability vector has different posterior distributions by treatment arm.
Each prior weight can be interpreted as the number of patients that have observed that 30-day response value.
We assume a-priori that each response value is equally likely with a total of \(6\) patients worth of data.
For the prior hyperparameters on the transition probability matrices we choose the identity matrix. This corresponds to assuming that there is one subject that has made the transition from each of the states to the same state, and that there are no subjects that changed state. Each early state vector carries a total weight of 1 subject, so this can be considered a weakly-informative prior.
# Priors for the longitudinal model
long_prior = list(
beta = rep(1, K), # priors for p30
alpha = diag(rep(1, K)) # priors for p90 | p30 and p180|p90
)We can visualize the priors with a heatmap of the possible transitions:
plt_prior = plot_Ptrans(list(p30 = long_prior$beta,
p90 = long_prior$alpha,
p180 = long_prior$alpha))
We can then run the fit_model() function and extract the predicted probabilities of success.
fit_pp = fit_model(df = data, long_prior = long_prior, n_max = 500,
key_trt = key_trt, iters = 1000, seed = model_seed)================================================================================
Imputation complete in 2.1 minutescat("The probability of success at the current sample size is", round(mean(fit_pp$PPn), 2), "\n")The probability of success at the current sample size is 1 cat("The probability of success at the maxiumum sample size is", round(mean(fit_pp$PPmax), 2), "\n")The probability of success at the maxiumum sample size is 0.98 The threshold to declare early success is generally based on \(PP_{n}\) and is chosen in the design phase to control the overall type I error. The threshold to declare futility is generally based on \(PP_{max}\) and is chosen in the design phase to minimize the exposure of patient to non-promising therapies.
In this case, we can see that the predictive probabilities of success if we stop now and if we continue to the maximum sample size are really high, despite the fact that the current p-value is not significant. This is due to the fact that the accruing information is expected to confirm the current positive results, reducing the uncertainty around the model estimates and yielding a significant p-value.
8 Bonus: Sparsity Inducing Prior for the Cutpoints
When the ordinal response categories are sparse (i.e., no observations are available for some categories), the model used is not likelihood identifiable. This is very common when the number of categories is large. In such scenario, using a uniform improper prior for the latent cutpoints would yield poor mixing and model degeneracies. Ideally, we would like to use domain expertise to inform a prior for the internal cutpoints. However, defining a prior on the abstract latent space where the internal cutpoints are defined can be difficult.
Say that we want to use a Dirichlet prior directly on the ordinal class probabilities. In order to construct the induced prior density function on the internal cutpoints we need to be careful about the Jacobian determinant of the transformation. In particular, the transformation from cutpoints to probabilities is \[ \begin{aligned} &p_{1} = 1 - \Pi(-c_{1}) \\ &p_{2} = \Pi(-c_{1}) - \Pi(-c_{2}) \\ &\vdots \\ &p_{K - 1} = \Pi(-c_{K - 2}) - \Pi(-c_{K - 1}) \\ &p_{K} = \Pi(-c_{K-1}) \end{aligned} \] with the constraint \(S = \sum_{k=1}^{K} p_{k} = 1\). The Jacobian corresponding to this transformation is \[ J = \left( \begin{array}{ccccc} 1 & \pi(-c_{1}) & & & \\ 1 & -\pi(-c_{1}) & \pi(-c_{2}) & & \\ \vdots & & \ddots & \ddots & \\ 1 & & & -\pi(-c_{k-2}) & \pi(-c_{k-1}) \\ 1 & & & & -\pi(-c_{k-1}) \end{array} \right) \]
Thus, the probability distribution on the cutpoints implied by the Dirichlet distribution on the probabilities is \[p(c_{1}, \dots, c_{K-1}) = p(p_{1}, \dots, p_{K} \mid \alpha) |J|.\]
The STAN model corresponding to this new prior is the following:
writeLines(readLines("./ord_logistic_dirichlet.stan"))functions {
real induced_dirichlet_lpdf(vector c, vector alpha) {
int D = num_elements(c) + 1;
vector[D - 1] sigma = inv_logit(- c);
vector[D] p;
matrix[D, D] J = rep_matrix(0, D, D);
// Induced ordinal probabilities
p[1] = 1 - sigma[1];
for (k in 2:(D - 1)) {
p[k] = sigma[k - 1] - sigma[k];
}
p[D] = sigma[D - 1];
// Baseline column of Jacobian
for (k in 1:D) {
J[k, 1] = 1;
}
// Diagonal entries of Jacobian
for (k in 2:D) {
J[k, k] = - sigma[k - 1] * (1 - sigma[k - 1]);
J[k - 1, k] = sigma[k - 1] * (1 - sigma[k - 1]);
}
return (dirichlet_lpdf(p | alpha) + log_determinant(J));
}
}
data {
int<lower=1> N; // Number of observations
int<lower=1> D; // Max of ordinal categories
int<lower=0, upper=1> X[N];
int<lower=1, upper=D> y[N]; // Observed ordinals {0, ..., K}
}
parameters {
real theta; // Latent effect
ordered[D - 1] c; // (Internal) cut points
}
model {
// Prior model
theta ~ normal(0, 2); // prior on the OR parameter
c ~ induced_dirichlet(rep_vector(1, D));
// Observational model
for (i in 1:N) {
y[i] ~ ordered_logistic(0.0, c + X[i] * theta);
}
}9 Conclusions
We discussed here the potential use of ordinal regression when dealing with ordinal outcomes. An alternative approach involves weighting the ordinal outcomes by means of a utility function that reflects the clinical meaning of the response values, and using a model for continuous responses.
Original Computing Environment
writeLines(readLines(file.path(Sys.getenv("HOME"), ".R/Makevars")))
CXX14FLAGS += -O3 -mtune=native -arch x86_64 -ftemplate-depth-256sessionInfo()R version 4.1.2 (2021-11-01)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] kableExtra_1.3.4 knitr_1.39 MCMCpack_1.6-3
[4] forcats_0.5.1 egg_0.4.5 gridExtra_2.3
[7] coda_0.19-4 bayesplot_1.9.0 rstan_2.21.5
[10] StanHeaders_2.21.0-7 RColorBrewer_1.1-3 latex2exp_0.9.4
[13] purrr_0.3.4 stringr_1.4.0 tidyr_1.2.0
[16] dplyr_1.0.9 pammtools_0.5.8 ggalluvial_0.12.3
[19] ggplot2_3.3.6 MASS_7.3-57
loaded via a namespace (and not attached):
[1] colorspace_2.0-3 ellipsis_0.3.2 ggridges_0.5.3
[4] rstudioapi_0.13 listenv_0.8.0 farver_2.1.0
[7] MatrixModels_0.5-0 prodlim_2019.11.13 fansi_1.0.3
[10] mvtnorm_1.1-3 xml2_1.3.3 codetools_0.2-18
[13] splines_4.1.2 Formula_1.2-4 jsonlite_1.8.0
[16] mcmc_0.9-7 compiler_4.1.2 httr_1.4.3
[19] backports_1.4.1 assertthat_0.2.1 Matrix_1.4-1
[22] fastmap_1.1.0 lazyeval_0.2.2 cli_3.3.0
[25] htmltools_0.5.2 quantreg_5.93 prettyunits_1.1.1
[28] tools_4.1.2 gtable_0.3.0 glue_1.6.2
[31] reshape2_1.4.4 posterior_1.2.2 Rcpp_1.0.8.3
[34] jquerylib_0.1.4 vctrs_0.4.1 svglite_2.1.0
[37] nlme_3.1-157 iterators_1.0.14 tensorA_0.36.2
[40] xfun_0.31 globals_0.15.0 ps_1.7.0
[43] rvest_1.0.2 lifecycle_1.0.1 future_1.26.1
[46] scales_1.2.0 inline_0.3.19 SparseM_1.81
[49] yaml_2.3.5 loo_2.5.1 sass_0.4.1
[52] stringi_1.7.6 highr_0.9 foreach_1.5.2
[55] checkmate_2.1.0 pkgbuild_1.3.1 lava_1.6.10
[58] rlang_1.0.2 pkgconfig_2.0.3 systemfonts_1.0.4
[61] matrixStats_0.62.0 distributional_0.3.0 evaluate_0.15
[64] lattice_0.20-45 labeling_0.4.2 processx_3.6.0
[67] tidyselect_1.1.2 parallelly_1.32.0 plyr_1.8.7
[70] magrittr_2.0.3 bookdown_0.26 R6_2.5.1
[73] generics_0.1.2 DBI_1.1.2 pillar_1.7.0
[76] withr_2.5.0 mgcv_1.8-40 abind_1.4-5
[79] survival_3.3-1 tibble_3.1.7 future.apply_1.9.0
[82] crayon_1.5.1 utf8_1.2.2 rmarkdown_2.14
[85] timereg_2.0.2 grid_4.1.2 callr_3.7.0
[88] digest_0.6.29 webshot_0.5.3 numDeriv_2016.8-1.1
[91] pec_2022.05.04 RcppParallel_5.1.5 stats4_4.1.2
[94] munsell_0.5.0 viridisLite_0.4.0 bslib_0.3.1