Hi there
I have been trying to use the original Stan code (provided by Harrell and Goodrich) for an unconstrained partial proportional odds model to analyse some simulated data mimicking a two-arm trial (with clusters removed):
// pointwise log-likelihood contributions
vector pw_log_lik(vector alpha, vector beta, matrix tau,
row_vector[] X, row_vector[] Z, int[] y) {
int N = size(X);
vector[N] out;
int k = max(y); // assumes all possible categories are observed
int j;
real cj;
real cj1;
for (n in 1:N) {
real eta = X[n] * beta;
j = y[n];
if (j == 1) cj = -( alpha[1] + eta );
else if (j == 2) cj = alpha[1] + eta;
else cj = alpha[j - 1] + eta + Z[n] * tau[ , j - 2];
if(j > 1 && j < k) cj1 = alpha[j] + eta + Z[n] * tau[ , j - 1];
if (j == 1 || j == k) out[n] = log_inv_logit(cj);
// else out[n] = log(1./(1. + exp(-cj)) - 1./(1. + exp(-cj1)));
else out[n] = log_diff_exp(-log1p_exp(- cj),
-log1p_exp(- cj1));
// else out[n] = log(-log1p_exp(-cj) + log1p_exp(-cj1));
}
return out;
}
// Pr(y == j)
matrix Pr(vector alpha, vector beta, matrix tau,
row_vector[] X, row_vector[] Z, int[] y) {
int N = size(X);
int k = max(y); // assumes all possible categories are observed
matrix[N, k] out;
for(n in 1:N) {
real eta = X[n] * beta ;
for(j in 1 : k) {
real cj;
real cj1;
if (j == 1) cj = -( alpha[1] + eta );
else if (j == 2) cj = alpha[1] + eta;
else cj = alpha[j - 1] + eta + Z[n] * tau[ , j - 2];
if(j > 1 && j < k) cj1 = alpha[j] + eta + Z[n] * tau[ , j - 1];
if (j == 1 || j == k) out[n, j] = log_inv_logit(cj);
//else out[n, j] = log(1./(1. + exp(-cj)) - 1./(1. + exp(-cj1)));
else out[n, j] = log_diff_exp(-log1p_exp(-cj),
-log1p_exp(-cj1));
// else out[n, j] = log(-log1p_exp(-cj) + log1p_exp(-cj1));
}
}
return exp(out);
}
}
data {
int<lower = 1> N; // number of observations
int<lower = 1> p; // number of predictors
int<lower = 1> q; // number of non-PO predictors in Z
matrix[N, p] X; // matrix of CENTERED predictors
matrix[N, q] Z; // matrix of CENTERED PPO predictors
int<lower = 2> k; // number of outcome categories
int<lower = 1, upper = k> y[N]; // outcome on 1 ... k
// prior standard deviations
real<lower = 0> sds;
real<lower = 0> sdsppo;
real<lower = 0> conc;
}
transformed data {
row_vector[p] Xr[N];
row_vector[q] Zr[N];
for (n in 1:N) Xr[n] = X[n, ];
for (n in 1:N) Zr[n] = Z[n, ];
}
parameters {
vector[p] beta; // coefficients on X
matrix[q, k - 2] tau; // coefficients on Z
simplex[k] pi; // category probabilities for a person w/ average predictors
}
transformed parameters {
vector[k - 1] alpha; // intercepts
vector[N] log_lik; // log-likelihood pieces
for (j in 2:k) alpha[j - 1] = logit(sum(pi[j:k])); // predictors are CENTERED
log_lik = pw_log_lik(alpha, beta, tau, Xr, Zr, y);
}
model {
target += log_lik;
target += dirichlet_lpdf(pi | rep_vector(conc, k));
target += normal_lpdf(beta | 0, sds);
for (j in 1:(k - 2)) target += normal_lpdf(tau[ , j] | 0, sdsppo);
}
Below is my code to generate the Stan inputs:
# Function to make stan inputs.
makeStanData <- function(dat)
{
z <- data.frame("y" = as.numeric(dat[["y"]]),
"X" = as.numeric(dat[["a"]]))
z <- z[["X"]] == 1 & z[["y"]] > 2
stanData <- list(N = nrow(dat),
y = as.numeric(dat[["y"]]),
X = as.matrix(dat[["a"]]),
Z = as.matrix(as.numeric(z)),
p = 1,
q = 1,
k = length(states),
p_par = rep(1, length(states)),
sds = 2.5,
sdsppo = 2.5,
conc = 1)
}
The model estimates for \tau and \beta do not match with what I expect (see below). For example, when generating the data assuming a common odds ratio of 1 (log-odds ratio of 0), the estimates are \beta = -1.34 and \tau[1,1] = 8.13 for a three category ordinal outcome, when they should be close to 0. How can I ensure that I obtain the correct/similar estimates using the original Stan model code? Is it to do with the fact that the predictors are centred, and I have to recover the estimates somehow?
fit$summary()
# A tibble: 4,009 × 10
variable mean median sd mad q5 q95 rhat
<chr> <num> <num> <num> <num> <num> <num> <num>
1 lp__ -1571. -1570. 1.44 1.23 -1574. -1569. 1.00
2 beta[1] -1.34 -1.34 0.106 0.105 -1.51 -1.16 1.00
3 tau[1,1] 8.13 8.02 0.936 0.883 6.81 9.82 1.00
4 pi[1] 0.140 0.140 0.0103 0.0102 0.124 0.158 1.00
5 pi[2] 0.536 0.536 0.0133 0.0131 0.515 0.558 1.00
6 pi[3] 0.323 0.323 0.0150 0.0151 0.299 0.348 1.00
7 alpha[1] 1.81 1.81 0.0851 0.0851 1.68 1.95 1.00
8 alpha[2] -0.740 -0.739 0.0687 0.0689 -0.854 -0.628 1.00
When compared to the model output using blrm with the exact same data, the output appears correct (i.e. blrm(singleDat$y~x, ppo=~x)) where both estimated cumulative log-odds ratios are close to 0 as expected (also below).
Mode Beta Mean Beta Median Beta S.E. Lower Upper Pr(Beta>0) Symmetry
y>=2 1.3615 1.3623 1.3631 0.0796 1.2047 1.5125 1.0000 1.02
y>=3 -0.5151 -0.5160 -0.5159 0.0639 -0.6389 -0.3898 0.0000 1.06
x 0.0950 0.0945 0.0972 0.1152 -0.1355 0.3163 0.7980 1.00
x:y>=3 -0.0314 -0.0305 -0.0305 0.1176 -0.2565 0.1993 0.4030 1.03
When I also run ppomod$rstan (where ppomod is a blrm object), I notice that the estimates for \beta and \tau are quite different to the original output. What is the difference between these estimates using $rstan?
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
beta[1] 2.11 0.04 2.58 -2.95 0.37 2.17 3.80 7.18 3320 1
tau[1,1] -0.68 0.05 2.63 -5.71 -2.43 -0.68 1.12 4.48 3301 1
alpha[1] 1.41 0.00 0.06 1.30 1.37 1.41 1.45 1.52 3808 1
alpha[2] -0.48 0.00 0.05 -0.57 -0.52 -0.48 -0.45 -0.39 3995 1
Here is a link to the Github repository below if this is helpful: the files simulateTrial and getProbs are both functions, and twoArmSim_example is the simulation/analysis file. Link to the repository is here: GitHub - chrisselman/ordpo
Ideally it would be great to use the original Stan code to generate things like posterior predictive distributions, as well as the posterior for each cumulative odds ratio which blrm does not currently have (I don’t think anyway - happy to be proven wrong!).
Any guidance would be much appreciated. :)