Ordinal Logistic Regression - Strange Results

I've been having struggling with the analysis of ordinal categorical data that I have (against 3 binary categorical explanatory variables). The problem is that for one of my data sets,running ordinal logistic regression in SPSS gives an unusual result, which raises questions about the validity of this statistical testing approach to my data: Full SPSS syntax and output can be accessed here

The question arises: why is no statistical result given for the non-redundant level of the age variable, StErr = 0, confidence interval does not exist, logit Estimate so high…? Exponentiation of logit gives some ridiculously high odds ratio… It’s like a perfect result with 100% certainty? Is this statistical test and outcome valid? Can someone explain what is happening in layman's terms? Should I be using a different statistical test than ordinal logistic regression? And please suggest how I should analyse and report it for a scientific publication.

The raw data can be viewed for further insight: I have also made available the SPSS data file in case you wanted to run the test for yourself:

Access here: SPSS .sav file

I have also attached a visualisation of the data in the form of stacked columns, to help you to quickly familiarise yourself with the character of my data. • can you explain what your variables mean and what you're trying to model here? Collagen level as a function of age, treatment, and genotype? What do the subscript 2 for these variables mean? Do they refer to age? You also have some very small cell sizes and if I understand right, empty cells (e.g. cells with 0 frequencies). – Marquis de Carabas Apr 6 '16 at 8:50
• Correct - collagen level as a function of age, treatment, and genotype. On that particular occasion that I generated that particular output, the subscript "2" for the explanatory variables is of no consequence - I just had to recode the original variables because I want to force SPSS to use the reference group that I wanted. The new resulting recoded variable is labelled "2". In the .sav data file I provided, there is no such "2" distinction - I just recoded into the original variable. Correct - small numbers, zero occurrences, and some purposely empty cells by design. – ptrcao Apr 6 '16 at 9:06
• Age does not seem to be modeled correctly. Age must be a continuous variable to give a correct statistical analysis otherwise you have unaccounted for heterogeneity in age. – Frank Harrell Apr 8 '16 at 13:02
• @FrankHarrell unaccounted heterogeneity? Could you please explain what you mean by this and how it arises with a binary age variable here? (Young = 3-5 months old, Old = 18-23 months old (mice).) – ptrcao Apr 9 '16 at 1:51
• I didn't realize that the ages were controlled to that extent. Unexplained heterogeneity in responses due to age would be present if there is any difference in response tendencies for 3m vs 5m old or 18m vs 23m old mice. If these age ranges made a difference, age would better be modeled as a continuous variable. – Frank Harrell Apr 9 '16 at 12:33

The problem here is what is called separation. For one level of your predictor (age) you have zero occurrences of one level of your outcome (collagen). The program is driving its estimate of the coefficient upwards to infinity and in this case stopping at an arbitrary value which as you point out when exponentiated gives you a very large odds ratio. What you do about it is really a substantive question to do with the underlying science rather than a statistical one.

You ask what separation means. In logistic regression it means that for some value of the predictor (or some linear combination of the predictors) only one value of the outcome occurs. You have a variant of this because you are using an ordinal model rather than logistic regression. What you do is a matter of your underlying science. Finding a perfect predictor may be really desirable for you or you may feel that this is just a chance finding caused by your sample size. There are models for logistic regression (exact logistic regression, Firth's bias corrected method) but I do not know whether they are available for ordinal models.

Investigating this further with two R packages (I do not use SPSS) reveals that they agree on there being a problem but give slightly different output. The function polr from MASS basically agrees with the output in the original question and inspecting the variance-covariance matrix of the coefficients shows that the estimated correlation between age and one of the intercepts is close to 1 in absolute value.

> fit.polr  summary(fit.polr)
Call:
polr(formula = coll ~ age + treat + gene, data = res, Hess = TRUE)

Coefficients:
Value Std. Error  t value
ageyoung -18.1584     0.2874 -63.1855
treatNMN  -0.4546     0.7457  -0.6096
geneig     0.4464     0.5823   0.7667

Intercepts:
Value    Std. Error t value
0|1 -17.4775   0.2871   -60.8811
1|2   0.7467   0.5709     1.3080

Residual Deviance: 68.77266
AIC: 78.77266
> cov2cor(vcov(fit.polr))
ageyoung    treatNMN     geneig         0|1        1|2
ageyoung  1.00000000 -0.06065359 -0.5308460 -0.99998399 -0.3224864
treatNMN -0.06065359  1.00000000  0.1169754  0.06384212  0.5512079
geneig   -0.53084598  0.11697536  1.0000000  0.53346172  0.6110366
0|1      -0.99998399  0.06384212  0.5334617  1.00000000  0.3278369
1|2      -0.32248636  0.55120792  0.6110366  0.32783686  1.0000000


The function clm from the ordinal package rejects the model on the grounds that the Hessian is singular. Using @FrankHarrell suggestion to investigate LRT we find that age is indeed significant.

> fit.clm  summary(fit.clm)
formula: coll ~ age + treat + gene
data:    res

logit flexible  54   -34.39 78.77 21(0) 5.63e-09 4.5e+09

Coefficients:
Estimate Std. Error z value Pr(>|z|)
ageyoung -23.3992         NA      NA       NA
treatNMN  -0.4546         NA      NA       NA
geneig     0.4465         NA      NA       NA

Threshold coefficients:
Estimate Std. Error z value
0|1 -22.7184         NA      NA
1|2   0.7467         NA      NA
> drop1(fit.clm, test = "Chisq")
Single term deletions

Model:
coll ~ age + treat + gene
Df     AIC    LRT  Pr(>Chi)
78.773
age     1 109.773 33.000 9.215e-09 ***
treat   1  77.151  0.378    0.5387
gene    1  77.365  0.592    0.4417
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Note the a LRT of age is possible here.

My feeling is that there is insufficient information here in the data-set to fit the model with any confidence. If we had a single observation for age == old & coll == 0 then both functions would fit the model (and give identical results) so the problem here is the zeroes for both (age == old & coll == 0) & (age == young & coll == 2).

• Thanks for articulating the problem clearly. So how should I analyse and report it for a scientific publication? – ptrcao Apr 7 '16 at 18:18
• And in what sense is the the term "separation" meant? Separation of what...? – ptrcao Apr 8 '16 at 5:16
• It's harder to have separation with ordinal models because of the order restriction. It would be good to double check this. If that is the problem the simplest fix is not to use Wald tests but to rely on likelihood ratio tests. – Frank Harrell Apr 8 '16 at 13:01
• So, just so I can understand, what is being separated from what...? – ptrcao Apr 8 '16 at 20:30
• @ptrcao see stats.stackexchange.com/questions/11109/… – AdamO Jan 4 '18 at 21:58