Ordinal regression with stan

Introduction

We extend the previous example on modelling the number of days where patients did not take their medication (grouped into bins):

less than 7 days
7-13 days
14-20 days
21-27 days
28+ days

accepting that the actual number of days that someone did not take medication could be 1.63, or 22.7 etc.

We introduce a continuous coviariate and also account for group-level heterogeneity.

Data

Generate some data for an ordinal outcome adheres to proportional odds and that aligns with the parameterisation used by stan. Initialise some reference cutpoints from which there will be deviations relating to dose assigned and group membership.

K <- 5
# probability of being in each category
p0 <- c(0.10, 0.30, 0.2, 0.35, 0.05)
p0[K] = 1 - sum(p0[1:(K-1)])
# cumulative probabilities
pc0 <- cumsum(p0)
# cumulative logits - these are the cutpoints for the control group
l0 <- qlogis(pc0[-K])
# shoud be the same 
rbind(p0, round(c(1, plogis(-l0[-K])) - c(plogis(-l0[-K]), 0), 3))

##    [,1] [,2] [,3] [,4] [,5]
## p0  0.1  0.3  0.2 0.35 0.05
##     0.1  0.3  0.2 0.35 0.05

Now assume a continous covariate and three groups introducing heterogeneity both in the intercepts and the rate of change attributable to the dose.

In other words, accommodate the possibility that there exists group-specific probabilities of each category and group specific changes that are associated with the dosage value. This will make more sense once we are looking at the model outputs.

set.seed(24099)
# 500 observations
N <- 1500
# group (1:3) level variation (random intercept)
u0_j <- rnorm(N, 0, 1)
u1_j <- rnorm(N, 0, 1)

grp_labs <- c("Ctl", "Trt1", "Trt2")
dose_options <- seq(0,1,by=0.1)

# log cumulative odds ratio for a unit increase in the covariate
b <- -1.1


get_y_obs <- function(x, j, b){
  # cumulative logits for treatment group
  l1 <- (b+u1_j[j])*x + u0_j[j] - l0
  # l0 is cutpoints, so transform l1 to get cumulative probs
  pc1 <- plogis(l1)
  p1 <- c(1, pc1) - c(pc1, 0)
  sample(1:K, 1, replace = F, p1)
}
get_y_obs <- Vectorize(get_y_obs)

# For each observation there is a dose (x) and they belong to one of
# three groups
d <- data.table(
  # x is distributed uniformly-ish
  # in theory, these could be any value across the range 
  # but here I just discretise to the dose_options
  x = dose_options[sample(1:length(dose_options), size = N, replace = T)],
  grp = sample(1:3, size = N, replace = T)
)
d[, y := get_y_obs(x, grp, b)]

# Calculate the true probabilities so that there is something
# with which to make a comparison.
get_prob <- function(x, j, k, b){
  # cumulative logits for treatment group
  l1 <- (b+u1_j[j])*x + u0_j[j] - l0
  # l0 is cutpoints, so transform l1 to get cumulative probs
  pc1 <- plogis(l1)
  p1 <- c(1, pc1) - c(pc1, 0)
  p1[k]
}
d_tru <- CJ(dose = dose_options, j = 1:3, k = 1:5)
d_tru$p_tru <- unlist(lapply(1:nrow(d_tru), function(i){
  get_prob(d_tru$dose[i], d_tru$j[i], d_tru$k[i], b)
}))
d_tru[, grp := grp_labs[j]]
d_tru[, q_tru := qlogis(p_tru)]

Model

The following model is a stan implementation. This model only includes the group level intercepts.

m1 <- cmdstanr::cmdstan_model(paste0(path_pref, "/stan/cumulative_logit4.stan"))
m1

## data{
##   int N;
##   int K;
##   int U; // num groups
##   array[N] int y;
##   vector[N] x;
##   vector[11] xc;
##   // group would be dp, np0, np1
##   array[N] int grp;
##   int prior_only;
##   real pri_b_s;
##   real pri_c_s;
##   real pri_u_s_r;
## }
## parameters {
##   real b;
##   ordered[K - 1] c;
##   vector[U] u0;
##   real<lower=0> u0_s;
## }
## model {
##   target += normal_lpdf(b | 0, pri_b_s);
##   target += normal_lpdf(c | 0, pri_c_s);
##   target += normal_lpdf(u0 | 0, 1);
##   target += exponential_lpdf(u0_s | pri_u_s_r);
##   if(!prior_only){
##     for (i in 1:N) {
##       target += ordered_logistic_lpmf(y[i] | x[i]*(b) + u0[grp[i]]*u0_s, c);
##     }
##   }
## }
## generated quantities{
##   array[K,11,U] real p;
##   
##   // want the predictions in each group. 
##   // groups (dp, np0, np1) 
##   for(i in 1:U){
##     // increments along the dose range
##     for(j in 1:11){
##       // categories
##       for(k in 1:K){
##         if(k == 1){
##           p[k,j,i] = 1 - inv_logit(xc[j]*(b) + u0[i]*u0_s - c[1]);
##         } else if (k == K) {
##           p[k,j,i] = inv_logit(xc[j]*(b) + u0[i]*u0_s - c[K-1]);
##         } else {
##           p[k,j,i] = inv_logit(xc[j]*(b) + u0[i]*u0_s - c[k-1]) - 
##             inv_logit(xc[j]*(b) + u0[i]*u0_s - c[k]);
##         }
##       }
##     }
##   } 
##   
## 
## }

Fit the model.

lsd <- list(
  N = nrow(d), K = K, U = 3,
  y = d$y, x = d$x, grp = d$grp,
  # predict effects for evenly spaced doses across the possible range
  xc = dose_options,
  prior_only = 0, pri_b_s = 1, pri_c_s = 1.8, pri_u_s_r = 1
)

f1 <- m1$sample(lsd,
    iter_warmup = 1000,
    iter_sampling = 2000,
    parallel_chains = 2,
    chains = 2,
    refresh = 500, show_exceptions = F)

## Running MCMC with 2 parallel chains...
## 
## Chain 1 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 2 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 2 Iteration:  500 / 3000 [ 16%]  (Warmup) 
## Chain 1 Iteration:  500 / 3000 [ 16%]  (Warmup) 
## Chain 1 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 1 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 2 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 2 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 1 Iteration: 1500 / 3000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 1500 / 3000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 2 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 1 Iteration: 2500 / 3000 [ 83%]  (Sampling) 
## Chain 2 Iteration: 2500 / 3000 [ 83%]  (Sampling) 
## Chain 1 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 1 finished in 140.4 seconds.
## Chain 2 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 2 finished in 155.4 seconds.
## 
## Both chains finished successfully.
## Mean chain execution time: 147.9 seconds.
## Total execution time: 155.5 seconds.

## Warning: 2 of 4000 (0.0%) transitions ended with a divergence.
## See https://mc-stan.org/misc/warnings for details.

f1$summary(variables = c("b", "u0_s", "u0"))

## # A tibble: 5 × 10
##   variable   mean median    sd   mad     q5    q95  rhat ess_bulk ess_tail
##   <chr>     <num>  <num> <num> <num>  <num>  <num> <num>    <num>    <num>
## 1 b        -1.01  -1.01  0.153 0.153 -1.26  -0.764  1.00    4329.    2431.
## 2 u0_s      1.34   1.20  0.563 0.475  0.679  2.46   1.00    1561.    1970.
## 3 u0[1]    -1.04  -1.00  0.543 0.550 -1.99  -0.217  1.00    1456.    2009.
## 4 u0[2]     1.09   1.01  0.663 0.663  0.136  2.29   1.00    1314.    1767.
## 5 u0[3]    -0.457 -0.450 0.456 0.455 -1.23   0.276  1.00    1423.    1878.

Within the model I use the generated quantities block to compute an estimate of the probability of observing each category at each dose within each group. Ordinarily, the group variable would be treated as a fixed effect, but in this example I am using what is commonly referred to as random effect.

post_p <- data.table(f1$draws(variables = c("p"), format = 'matrix'))
post_p <- melt(post_p, measure.vars = names(post_p))
# categories
post_p[, k := gsub("p[", "", variable, fixed = T)]
post_p[, k := gsub(",.*$", "", k)]
# dosage
post_p[, xi := gsub("^p\\[.,", "", variable)]
post_p[, xi := gsub(",.]", "", xi)]
post_p[, xi := factor(xi, levels = paste0(1:11))]
# group
post_p[, trt := gsub("^p\\[.,.*,", "", variable)]
post_p[, trt := gsub("]", "", trt, fixed = T)]
post_p[trt == 1, grp := "Ctl"]
post_p[trt == 2, grp := "Trt1"]
post_p[trt == 3, grp := "Trt2"]

post_p[, dose := lsd$xc[xi]]

d_tbl <- post_p[, .(p_mu = mean(value)), keyby = .(grp, k, dose)]
d_tbl <- melt(d_tbl, id.vars = c("grp", "k", "dose"))

d_tbl <- dcast(
  d_tbl, dose + grp + variable ~ k, value.var = "value"
)
# compare to the truth
knitr::kable(d_tbl, format = "simple", digits = 3)

dose	grp	variable	1	2	3	4	5
0.0	Ctl	p_mu	0.319	0.401	0.141	0.124	0.015
0.0	Trt1	p_mu	0.039	0.144	0.167	0.500	0.149
0.0	Trt2	p_mu	0.194	0.375	0.191	0.211	0.029
0.1	Ctl	p_mu	0.341	0.399	0.132	0.114	0.014
0.1	Trt1	p_mu	0.043	0.156	0.175	0.489	0.137
0.1	Trt2	p_mu	0.210	0.383	0.184	0.195	0.026
0.2	Ctl	p_mu	0.364	0.395	0.124	0.104	0.012
0.2	Trt1	p_mu	0.048	0.168	0.182	0.477	0.126
0.2	Trt2	p_mu	0.227	0.390	0.177	0.181	0.024
0.3	Ctl	p_mu	0.388	0.389	0.116	0.096	0.011
0.3	Trt1	p_mu	0.052	0.180	0.189	0.464	0.115
0.3	Trt2	p_mu	0.246	0.396	0.170	0.167	0.022
0.4	Ctl	p_mu	0.412	0.382	0.109	0.087	0.010
0.4	Trt1	p_mu	0.058	0.194	0.195	0.449	0.105
0.4	Trt2	p_mu	0.265	0.400	0.162	0.154	0.020
0.5	Ctl	p_mu	0.436	0.373	0.101	0.080	0.009
0.5	Trt1	p_mu	0.063	0.207	0.201	0.433	0.096
0.5	Trt2	p_mu	0.285	0.402	0.154	0.142	0.018
0.6	Ctl	p_mu	0.461	0.363	0.094	0.073	0.008
0.6	Trt1	p_mu	0.069	0.221	0.206	0.416	0.087
0.6	Trt2	p_mu	0.306	0.402	0.146	0.131	0.016
0.7	Ctl	p_mu	0.487	0.352	0.087	0.067	0.008
0.7	Trt1	p_mu	0.076	0.236	0.209	0.399	0.080
0.7	Trt2	p_mu	0.327	0.401	0.137	0.120	0.015
0.8	Ctl	p_mu	0.512	0.340	0.081	0.061	0.007
0.8	Trt1	p_mu	0.084	0.251	0.212	0.381	0.073
0.8	Trt2	p_mu	0.350	0.397	0.129	0.110	0.013
0.9	Ctl	p_mu	0.537	0.327	0.074	0.055	0.006
0.9	Trt1	p_mu	0.092	0.265	0.214	0.362	0.066
0.9	Trt2	p_mu	0.373	0.393	0.121	0.101	0.012
1.0	Ctl	p_mu	0.562	0.314	0.069	0.050	0.006
1.0	Trt1	p_mu	0.101	0.280	0.215	0.344	0.060
1.0	Trt2	p_mu	0.397	0.386	0.113	0.093	0.011

The following plots show the probability of observing each category in each group. I have discreticised the dose values into chunks of 0.1 for the pursposes of visualisation.

d_fig <- post_p[, .(p_mu = mean(value)), keyby = .(k, dose, grp)]
ggplot(d_fig, aes(x = dose, y = p_mu, fill = k)) +
  geom_bar(position="fill", stat="identity") +
  scale_x_continuous("dose", breaks = c(0, 0.5, 1)) +
  scale_fill_discrete("Category") +
  facet_wrap(~grp) +
  theme_bw() +
  theme(legend.position = "bottom")

And in the following plot, I show the posterior means for the log-odds of each category at each dose in each group. You can see how the dosage slopes are quite different for each group and the model fit is not very good. Also, I show the true log-odds from which the data was simulated as points.

d_fig <- post_p[, .(
  q_mu = mean(qlogis(value)),
  q_025 = quantile(qlogis(value), prob = 0.025),
  q_975 = quantile(qlogis(value), prob = 0.975)
  ), keyby = .(k, dose, grp)]

ggplot(d_fig, aes(x = dose, y = q_mu, col = grp)) +
  geom_ribbon(aes(ymin = q_025, ymax = q_975, 
                  fill = grp, col = NULL),
              alpha = 0.1) +
  geom_line() +
  geom_point(data = d_tru,
             aes(x = dose, y = q_tru, col = grp), 
             size = 0.4) +
  scale_x_continuous("Dose", breaks = seq(0,1,by=0.2)) +
  scale_y_continuous("Pr(Y = y)") +
  scale_color_discrete("") +
  scale_fill_discrete("") +
  facet_wrap(~paste0("Category ", k)) +
  theme_bw() +
  theme(legend.position = "bottom")

The next figure shows the same estimates on the probability scale (along with the true probabilities as points).

d_fig <- post_p[, .(
  p_mu = mean(value),
  p_025 = quantile(value, prob = 0.025),
  p_975 = quantile(value, prob = 0.975)
  ), keyby = .(k, dose, grp)]

d[, id_grp := .GRP, keyby = .(grp, x)]
d[, total := .N, keyby = id_grp]
d_tmp <- d[, .(p_obs = .N/total), keyby = .(grp, x, y)]
setnames(d_tmp, "x", "dose")
setnames(d_tmp, "y", "k")
d_tmp$k <- factor(d_tmp$k)
d_tmp$grp <- factor(d_tmp$grp, levels = 1:3, labels = c(
  "Ctl", "Trt1", "Trt2"
))

ggplot(d_fig, aes(x = dose, y = p_mu, col = grp)) +
  geom_ribbon(aes(ymin = p_025, ymax = p_975, 
                  fill = grp, col = NULL),
              alpha = 0.1) +
  geom_line() +
  geom_point(data = d_tru,
             aes(x = dose, y = p_tru, col = grp), 
             size = 0.4) +
  scale_x_continuous("Dose", breaks = seq(0,1,by=0.2)) +
  scale_y_continuous("Pr(Y = y)") +
  scale_color_discrete("") +
  scale_fill_discrete("") +
  facet_wrap(~paste0("Category ", k)) +
  theme_bw() +
  theme(legend.position = "bottom")

Now introduce the group level adjustment to the dose effect as well as the random intercept.

m2 <- cmdstanr::cmdstan_model(paste0(path_pref, "/stan/cumulative_logit5.stan"))
m2

## data{
##   int N;
##   int K;
##   int U; // num groups
##   array[N] int y;
##   vector[N] x;
##   vector[11] xc;
##   // group would be dp, np0, np1
##   array[N] int grp;
##   int prior_only;
##   real pri_b_s;
##   real pri_c_s;
##   real pri_u_s_r;
## }
## parameters {
##   real b;
##   ordered[K - 1] c;
##   vector[U] u0;
##   real<lower=0> u0_s;
##   vector[U] u1;
##   real<lower=0> u1_s;
## }
## model {
##   target += normal_lpdf(b | 0, pri_b_s);
##   target += normal_lpdf(c | 0, pri_c_s);
##   target += normal_lpdf(u0 | 0, 1);
##   target += exponential_lpdf(u0_s | pri_u_s_r);
##   target += normal_lpdf(u1 | 0, 1);
##   target += exponential_lpdf(u1_s | pri_u_s_r);
##   if(!prior_only){
##     for (i in 1:N) {
##       target += ordered_logistic_lpmf(y[i] | x[i]*(b + u1[grp[i]]*u1_s) + u0[grp[i]]*u0_s, c);
##     }
##   }
## }
## generated quantities{
##   array[K,11,U] real p;
##   
##   // want the predictions in each group. 
##   // groups (dp, np0, np1) 
##   for(i in 1:U){
##     // increments along the dose range
##     for(j in 1:11){
##       // categories
##       for(k in 1:K){
##         if(k == 1){
##           p[k,j,i] = 1 - inv_logit(xc[j]*(b+u1[i]*u1_s) + u0[i]*u0_s - c[1]);
##         } else if (k == K) {
##           p[k,j,i] = inv_logit(xc[j]*(b+u1[i]*u1_s) + u0[i]*u0_s - c[K-1]);
##         } else {
##           p[k,j,i] = inv_logit(xc[j]*(b+u1[i]*u1_s) + u0[i]*u0_s - c[k-1]) - 
##             inv_logit(xc[j]*(b+u1[i]*u1_s) + u0[i]*u0_s - c[k]);
##         }
##       }
##     }
##   } 
##   
## 
## }

Fit the random intercept and random slope model (assuming independence between the two).

lsd <- list(
  N = nrow(d), K = K, U = 3,
  y = d$y, x = d$x, grp = d$grp,
  # predict effects for evenly spaced doses across the possible range
  xc = dose_options,
  prior_only = 0, pri_b_s = 1, pri_c_s = 1.8, pri_u_s_r = 1
)

f2 <- m2$sample(lsd,
    iter_warmup = 1000,
    iter_sampling = 2000,
    parallel_chains = 2,
    chains = 2,
    refresh = 500, show_exceptions = F)

## Running MCMC with 2 parallel chains...
## 
## Chain 1 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 2 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 2 Iteration:  500 / 3000 [ 16%]  (Warmup) 
## Chain 1 Iteration:  500 / 3000 [ 16%]  (Warmup) 
## Chain 2 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 2 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 1 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 1 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 2 Iteration: 1500 / 3000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 1500 / 3000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 1 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 2 Iteration: 2500 / 3000 [ 83%]  (Sampling) 
## Chain 1 Iteration: 2500 / 3000 [ 83%]  (Sampling) 
## Chain 2 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 2 finished in 161.6 seconds.
## Chain 1 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 1 finished in 181.6 seconds.
## 
## Both chains finished successfully.
## Mean chain execution time: 171.6 seconds.
## Total execution time: 181.7 seconds.

## Warning: 2 of 4000 (0.0%) transitions ended with a divergence.
## See https://mc-stan.org/misc/warnings for details.

f2$summary(variables = c("b", "u0_s", "u0", "u1_s", "u1"))

## # A tibble: 9 × 10
##   variable   mean median    sd   mad     q5     q95  rhat ess_bulk ess_tail
##   <chr>     <num>  <num> <num> <num>  <num>   <num> <num>    <num>    <num>
## 1 b        -0.674 -0.718 0.642 0.594 -1.62   0.422   1.00    1747.    1993.
## 2 u0_s      0.945  0.817 0.492 0.377  0.415  1.89    1.00    1971.    1946.
## 3 u0[1]    -0.888 -0.850 0.583 0.591 -1.90   0.0150  1.00    1865.    2575.
## 4 u0[2]     0.907  0.853 0.694 0.698 -0.117  2.13    1.00    2191.    2579.
## 5 u0[3]    -0.431 -0.411 0.530 0.514 -1.32   0.426   1.00    1859.    2549.
## 6 u1_s      1.34   1.20  0.640 0.512  0.615  2.60    1.00    2005.    2664.
## 7 u1[1]    -1.13  -1.08  0.585 0.578 -2.13  -0.266   1.00    2281.    2501.
## 8 u1[2]     0.854  0.791 0.673 0.665 -0.128  2.03    1.00    1675.    2433.
## 9 u1[3]    -0.522 -0.500 0.513 0.491 -1.39   0.284   1.00    2020.    2431.

post_p <- data.table(f2$draws(variables = c("p"), format = 'matrix'))
post_p <- melt(post_p, measure.vars = names(post_p))
# categories
post_p[, k := gsub("p[", "", variable, fixed = T)]
post_p[, k := gsub(",.*$", "", k)]
# dosage
post_p[, xi := gsub("^p\\[.,", "", variable)]
post_p[, xi := gsub(",.]", "", xi)]
post_p[, xi := factor(xi, levels = paste0(1:11))]
# group
post_p[, trt := gsub("^p\\[.,.*,", "", variable)]
post_p[, trt := gsub("]", "", trt, fixed = T)]
post_p[trt == 1, grp := "Ctl"]
post_p[trt == 2, grp := "Trt1"]
post_p[trt == 3, grp := "Trt2"]

post_p[, dose := lsd$xc[xi]]

d_tbl <- post_p[, .(p_mu = mean(value)), keyby = .(grp, k, dose)]
d_tbl <- melt(d_tbl, id.vars = c("grp", "k", "dose"))

d_tbl <- dcast(
  d_tbl, dose + grp + variable ~ k, value.var = "value"
)
# compare to the truth
knitr::kable(d_tbl, format = "simple", digits = 3)

dose	grp	variable	1	2	3	4	5
0.0	Ctl	p_mu	0.222	0.397	0.177	0.179	0.024
0.0	Trt1	p_mu	0.067	0.223	0.206	0.416	0.089
0.0	Trt2	p_mu	0.167	0.366	0.200	0.233	0.034
0.1	Ctl	p_mu	0.259	0.407	0.162	0.153	0.019
0.1	Trt1	p_mu	0.065	0.219	0.205	0.420	0.091
0.1	Trt2	p_mu	0.186	0.380	0.193	0.211	0.029
0.2	Ctl	p_mu	0.299	0.410	0.146	0.129	0.016
0.2	Trt1	p_mu	0.064	0.216	0.204	0.424	0.093
0.2	Trt2	p_mu	0.207	0.392	0.184	0.191	0.026
0.3	Ctl	p_mu	0.343	0.406	0.129	0.108	0.013
0.3	Trt1	p_mu	0.062	0.213	0.203	0.428	0.094
0.3	Trt2	p_mu	0.230	0.401	0.174	0.173	0.023
0.4	Ctl	p_mu	0.390	0.395	0.113	0.091	0.011
0.4	Trt1	p_mu	0.061	0.209	0.202	0.431	0.096
0.4	Trt2	p_mu	0.254	0.407	0.164	0.155	0.020
0.5	Ctl	p_mu	0.439	0.378	0.098	0.076	0.009
0.5	Trt1	p_mu	0.059	0.206	0.201	0.435	0.098
0.5	Trt2	p_mu	0.280	0.410	0.153	0.139	0.017
0.6	Ctl	p_mu	0.489	0.356	0.084	0.063	0.007
0.6	Trt1	p_mu	0.058	0.203	0.200	0.439	0.100
0.6	Trt2	p_mu	0.308	0.410	0.142	0.125	0.015
0.7	Ctl	p_mu	0.540	0.330	0.072	0.052	0.006
0.7	Trt1	p_mu	0.057	0.200	0.198	0.442	0.103
0.7	Trt2	p_mu	0.337	0.407	0.131	0.111	0.013
0.8	Ctl	p_mu	0.589	0.302	0.061	0.043	0.005
0.8	Trt1	p_mu	0.056	0.197	0.197	0.446	0.105
0.8	Trt2	p_mu	0.367	0.401	0.121	0.099	0.012
0.9	Ctl	p_mu	0.637	0.272	0.051	0.036	0.004
0.9	Trt1	p_mu	0.055	0.194	0.195	0.449	0.107
0.9	Trt2	p_mu	0.399	0.392	0.110	0.088	0.010
1.0	Ctl	p_mu	0.682	0.242	0.043	0.030	0.003
1.0	Trt1	p_mu	0.054	0.191	0.194	0.452	0.110
1.0	Trt2	p_mu	0.431	0.381	0.101	0.079	0.009

Now, replot the posterior means on the log-odds and probability scales. The true probabilities are quite well recovered.

d_fig <- post_p[, .(
  q_mu = mean(qlogis(value)),
  q_025 = quantile(qlogis(value), prob = 0.025),
  q_975 = quantile(qlogis(value), prob = 0.975)
  ), keyby = .(k, dose, grp)]

ggplot(d_fig, aes(x = dose, y = q_mu, col = grp)) +
  geom_ribbon(aes(ymin = q_025, ymax = q_975, 
                  fill = grp, col = NULL),
              alpha = 0.1) +
  geom_line() +
  geom_point(data = d_tru,
             aes(x = dose, y = q_tru, col = grp), 
             size = 0.4) +
  scale_x_continuous("Dose", breaks = seq(0,1,by=0.2)) +
  scale_y_continuous("Pr(Y = y)") +
  scale_color_discrete("") +
  scale_fill_discrete("") +
  facet_wrap(~paste0("Category ", k)) +
  theme_bw() +
  theme(legend.position = "bottom")

d_fig <- post_p[, .(
  p_mu = mean(value),
  p_025 = quantile(value, prob = 0.025),
  p_975 = quantile(value, prob = 0.975)
  ), keyby = .(k, dose, grp)]

d[, id_grp := .GRP, keyby = .(grp, x)]
d[, total := .N, keyby = id_grp]
d_tmp <- d[, .(p_obs = .N/total), keyby = .(grp, x, y)]
setnames(d_tmp, "x", "dose")
setnames(d_tmp, "y", "k")
d_tmp$k <- factor(d_tmp$k)
d_tmp$grp <- factor(d_tmp$grp, levels = 1:3, labels = c(
  "Ctl", "Trt1", "Trt2"
))

ggplot(d_fig, aes(x = dose, y = p_mu, col = grp)) +
  geom_ribbon(aes(ymin = p_025, ymax = p_975, 
                  fill = grp, col = NULL),
              alpha = 0.1) +
  geom_line() +
  geom_point(data = d_tru,
             aes(x = dose, y = p_tru, col = grp), 
             size = 0.4) +
  scale_x_continuous("Dose", breaks = seq(0,1,by=0.2)) +
  scale_y_continuous("Pr(Y = y)") +
  scale_color_discrete("") +
  scale_fill_discrete("") +
  facet_wrap(~paste0("Category ", k)) +
  theme_bw() +
  theme(legend.position = "bottom")

Ordinal regression with stan

2023-08-24

Introduction

Data

Model