Different results using partial proportional odds model

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. :)