# Different results in unconstrained partial proportional odds model

**URL:** <https://discourse.datamethods.org/t/different-results-in-unconstrained-partial-proportional-odds-model/6798>\
**Category:** modeling strategy\
**Tags:** ordinal, bayes, rct\
**Created:** [May 22, 2023, 7:14am UTC](https://discourse.datamethods.org/t/different-results-in-unconstrained-partial-proportional-odds-model/6798 "2023-05-22T07:14:19Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![cselman](https://discourse.datamethods.org/letter_avatar_proxy/v4/letter/c/ccd318/32.png) [@cselman](https://discourse.datamethods.org/u/cselman)\
**Post date:** [May 22, 2023, 7:14am UTC](https://discourse.datamethods.org/t/different-results-in-unconstrained-partial-proportional-odds-model/6798/1 "2023-05-22T07:14:20Z")

</div>

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

```plaintext
  // 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:

```auto
# 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?

```auto
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).

```auto
       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?

```auto
             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](https://github.com/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. 🙂