Hi all,
Happy thanksgiving to those celebrating this weekend. I am doing some work extending/adapting @f2harrell ImpactPO approach to network meta-analysis or ordinal outcomes. As part of this workflow I am visualizing treatment effects on thresholds of an ordinal outcome + the OR from a PO model for each trial in the meta-analysis. I’ve generally thought of the OR from a PO model as a sort of average OR but have been surprised to find that in some studies the estimate of the PO OR can lie entirely outside of ORs from the individual thresholds. I’ve attached a 1 study reprex where this occurs below (ordinal risankizumab lor of 4.2 is greater than either pasi75 4.1 or pasi90 4.08). The pasi75/pasi90 outcomes are generated from the categorical version of the outcome, so I’m pretty sure it’s not just an issue with the data but would be grateful if anyone could help me realize a mistake a making. I haven’t been able to find anything talking about this on a quick scan of the literature so had a few questions to the more technically minded in our group:
- Is this a known phenomenon and if so is there an intuitive explanation? Is it related to non collapsibility at all?
- I’ve noticed it generally doesn’t matter when converting back to absolute scales in my case, but is there a good way to explain this to people who are already skeptical of using PO models? The deviations are usually small when it does happen, but for treatments with extreme treatment effects (lor > 7) it can start to look like very large differences on the OR scale (but then <0.1% differences on risk difference).
- @f2harrell I noticed in ImpactPO you use multinomial regression when relaxing PO on all variables at once. I noticed sometimes when fitting this same model with VGAM and parallel = FALSE for the treatment variable that I get a solution that converges but to a very bizarre solution where both the intercepts and treatment effects are not at all similar to running models on the individual thresholds separately. Do you know why this would be?
Example code and results
library(dplyr)
library(rms)
dat <- readxl::read_xlsx(here::here("ult1.xlsx"))
dat %>%
dplyr::group_by(trt) %>%
dplyr::summarise(p75 = mean(pasi75),
p90 = mean(pasi90),
n = dplyr::n())
#> # A tibble: 3 × 4
#> trt p75 p90 n
#> <fct> <dbl> <dbl> <int>
#> 1 Placebo 0.0980 0.0490 102
#> 2 Risankizumab 0.868 0.753 304
#> 3 Ustekinumab 0.7 0.420 100
p75 <- glm(pasi75 ~ trt, family = binomial, data = dat)
p90 <- glm(pasi90 ~ trt, family = binomial, data = dat)
ord <- lrm(out ~ trt, data = dat)
p75
#>
#> Call: glm(formula = pasi75 ~ trt, family = binomial, data = dat)
#>
#> Coefficients:
#> (Intercept) trtRisankizumab trtUstekinumab
#> -2.219 4.106 3.067
#>
#> Degrees of Freedom: 505 Total (i.e. Null); 503 Residual
#> Null Deviance: 634.5
#> Residual Deviance: 424.3 AIC: 430.3
p90
#>
#> Call: glm(formula = pasi90 ~ trt, family = binomial, data = dat)
#>
#> Coefficients:
#> (Intercept) trtRisankizumab trtUstekinumab
#> -2.965 4.082 2.642
#>
#> Degrees of Freedom: 505 Total (i.e. Null); 503 Residual
#> Null Deviance: 697.3
#> Residual Deviance: 515.6 AIC: 521.6
ord
#> Logistic Regression Model
#>
#> lrm(formula = out ~ trt, data = dat)
#>
#> Model Likelihood Discrimination Rank Discrim.
#> Ratio Test Indexes Indexes
#> Obs 506 LR chi2 234.27 R2 0.434 C 0.776
#> PASI < 75 162 d.f. 2 R2(2,506)0.368 Dxy 0.552
#> PASI 75 - 89 68 Pr(> chi2) <0.0001 R2(2,406.1)0.436 gamma 0.796
#> PASI 90-100 276 Brier 0.128 tau-a 0.322
#> max |deriv| 3e-12
#>
#> Coef S.E. Wald Z Pr(>|Z|)
#> y>=PASI 75 - 89 -2.2116 0.3327 -6.65 <0.0001
#> y>=PASI 90-100 -3.1152 0.3472 -8.97 <0.0001
#> trt=Risankizumab 4.2124 0.3679 11.45 <0.0001
#> trt=Ustekinumab 2.9048 0.3831 7.58 <0.0001