I am quite puzzled by the logistic regression results with three outcome categories (0,1,2); 0 is no feelings, 1 is slightly happy, 2 is extremely happy.

I tried both (1) logistic regression and ordered the outcome (2) using ordinal logistic regression through MASS::polr

The summary from (1) looks like this:

glm(formula = FeelingOutcome ~ Dosage + Age + factor(Sex.x) + 
    factor(Race.x) + TestPeriod, family = binomial, data = TestSet, 
    na.action = "na.exclude")

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.7767  -1.1058   0.6135   1.0088   2.0481  

                      Estimate Std. Error z value Pr(>|z|)  
(Intercept)          1.457e+01  1.455e+03   0.010   0.9920  
Dosage              -1.981e+00  8.145e-01  -2.433   0.0150 *
Age                  4.434e-02  2.494e-02   1.778   0.0755 .
factor(Sex.x)Male    6.504e-01  4.544e-01   1.431   0.1523  
factor(Race.x)Black -1.670e+01  1.455e+03  -0.011   0.9908  
factor(Race.x)White -1.513e+01  1.455e+03  -0.010   0.9917  
TestPeriod          -1.413e-04  1.319e-04  -1.072   0.2839  
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 146.34  on 105  degrees of freedom
Residual deviance: 130.71  on  99  degrees of freedom
  (2 observations deleted due to missingness)
AIC: 144.71

Number of Fisher Scoring iterations: 14

The statistics showed moderate significance of Dosage association with Feeling Outcome.

However, when I use ordinal logistic regression, the strange number showed:

  TestORdinalLR <- MASS::polr(FeelingOutcome ~ Dosage + 
                                     Age + 
                                     factor(Sex.x) +
                                     factor(Race.x) + 
                                   data=TestSet, Hess = TRUE, na.action = "na.exclude")
ctable <- coef(summary(TestORdinalLR))
p_OrdiDxSHCec <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
ctable <- cbind(ctable, "p value" = p_OrdiDxSHCec)
    > ctable
                                  Value   Std. Error      t value       p value
Dosage              -1.270935e+00 0.0948121344  -13.4047742  5.669557e-41
Age                  3.844785e-02 0.0138925651    2.7675126  5.648587e-03
factor(Sex.x)Male    1.799269e-01 0.3784498549    0.4754313  6.344796e-01
factor(Race.x)Black -1.670944e+01 0.3947162775  -42.3327766  0.000000e+00
factor(Race.x)White -1.536035e+01 0.4342894383  -35.3689268 5.131513e-274
TestPeriod          -4.578557e-05 0.0003773869   -0.1213226  9.034355e-01
0|1                 -1.432219e+01 0.0395987756 -361.6826279  0.000000e+00
1|2                 -1.365826e+01 0.1579986481  -86.4454177  0.000000e+00

I saw a whopping p-value change from 0.0150 to 5.669557e-41. Intuitively, I know I should use ordinal logistic regression, but from the results, the logistic regression seems more realistic?

Here's the data distribution (I flipped the x and y for visualization): enter image description here

  • 1
    $\begingroup$ DId you discretize the outcome somehow for your logistic regression? You didn't mention it. I would've expected it to error if you tried lgoistic regression with FeelingOutcome that has values outside of $[0,1]$. $\endgroup$
    – bdeonovic
    Commented Feb 18, 2022 at 19:59
  • 2
    $\begingroup$ A glm model with family=binomial will fit a logistic regression with two outcome levels only, not three outcome levels. $\endgroup$ Commented Feb 18, 2022 at 20:39
  • 1
    $\begingroup$ What is your reference level for factor(Race.x)? You show coefficients for both Black and White, so the reference level presumably has relatively few members in your data. You didn't get an error in binomial logistic regression because, as the help page says, "the response can also be specified as a factor (when the first level denotes failure and all others success)." Your outcome was presumably interpreted as a factor. $\endgroup$
    – EdM
    Commented Feb 18, 2022 at 20:58
  • 2
    $\begingroup$ I don't know why there isn't a warning. ... It appears that if there are more than two levels in the dependent variable, that glm(family=binomal) will use the first level as the first level, and all other levels as the second level. $\endgroup$ Commented Feb 18, 2022 at 21:01
  • 2
    $\begingroup$ I will read through some posts on that before I do anything else. For now, I think the question I asked is solved, it's pretty obvious to all of you that family=binomial in glm should only have two outcomes. Anything more than three should be multivariable or ordinal, and I will study a bit to see if I need more help. Thank you all for your time. $\endgroup$
    – Molly_K
    Commented Feb 18, 2022 at 21:14

1 Answer 1


As a quick check, ensure the variable is a ordered factor using str(TestSet). I do not know if this would be a problem for the polr function or not.

I think you may be using the analysis from here: https://stats.oarc.ucla.edu/r/dae/ordinal-logistic-regression/. If not then it is a good resource and explains things well. Their data set is bigger than yours, thus the manual p-value calculation you have done may be less satisfactory on your data set.

From the webpage: “One way to calculate a p-value in this case is by comparing the t-value against the standard normal distribution, like a z test. Of course this is only true with infinite degrees of freedom, but is reasonably approximated by large samples, becoming increasingly biased as sample size decreases.”.

The calculation you did for p-values matches theirs. However, your lower sample size means your p-values will be biased to be low.

Unless you have done so already checking the proportional odds assumption is worth doing. Per the webpage section proportional odds assumption: “One of the assumptions underlying ordinal logistic (and ordinal probit) regression is that the relationship between each pair of outcome groups is the same. In other words, ordinal logistic regression assumes that the coefficients that describe the relationship between, say, the lowest versus all higher categories of the response variable are the same as those that describe the relationship between the next lowest category and all higher categories, etc. This is called the proportional odds assumption or the parallel regression assumption. Because the relationship between all pairs of groups is the same, there is only one set of coefficients.”.

How to go about this is also detailed.

Hope this helps.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.