Ordinal logistic regression in R - p value interpretation I want to know if a certain protein can be a marker of a particular disease and I'm searching for a dose response.
I've created a toy example with three ordered groups: control, sick and very sick.
Although I purposely create a protein which doesn't follow a dose reponse, with clearly no overlap between groups, the ordinal logistic regression model outputs a significant coefficient for protein.

*

*Why?

*Which model could I use to only get true dose response, negative or positive?
Thank you.

library(ordinal)
library(ggplot2)

d=data.frame(
  group=rep(c('control','sick','very_sick'),times=c(205,205,205)),
  protein=c(rnorm(205, mean=10, sd=1), rnorm(205, mean=-5, sd=1), rnorm(205, mean=5, sd=1)))
 
d$group<- factor(d$group,levels = c('control','sick','very_sick'), ordered = TRUE)  
 

ggplot(d, aes(x=group, y=protein)) + 
  geom_boxplot()

model=clm(group ~ protein, data=d, link='logit')

summary(model)


## formula: group ~ protein
## data:    d
## 
##  link  threshold nobs logLik  AIC     niter max.grad cond.H 
##  logit flexible  615  -646.06 1298.12 4(0)  1.12e-11 1.0e+02
## 
## Coefficients:
##         Estimate Std. Error z value Pr(>|z|)    
## protein -0.08753    0.01152  -7.596 3.07e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Threshold coefficients:
##             Estimate Std. Error z value
## control|low -0.87181    0.09318  -9.356
## low|high     0.62542    0.08807   7.102



Why?

To extend on the comment by @sal-mangiafico
This is because an inadvertent negative trend was introduced.  this is easiest to see in this replication with graph below:
set.seed(11112022)

d=data.frame(
  group=rep(c('control','sick','very_sick'),times=c(205,205,205)),
  protein=c(rnorm(205, mean=10, sd=1), rnorm(205, mean=-5, sd=1), rnorm(205, mean=5, sd=1)))

d$group <- factor(d$group,levels = c('control','sick','very_sick'), ordered = TRUE) 

d$group2 <- as.numeric(d$group)

ggplot(d, aes(x=group, y=protein)) + 
  geom_boxplot() + stat_smooth(aes(x = group2), method = "lm")

Producing:

Because the outcome is treated as ordinal, the line fit has a slope given the pattern of means that trend downward given the implied ordering of the levels.
There is no problem with the model applied in this case.  The issue is, most likely, an unexpected result from the way the data were generated.
