I have been using ordinal logistic regression with the
ordinal package in R and the clm() function to complete an analysis of ordinal survey research data. In essence, I would simply like to discern the strength and direction of relationships between the predictor variables and the response variable (all of which are ordinal). Building a model to predict with is not the goal here.
Having already identified a set of 5 predictor variables from a possible 94, I have passed them on to an ordinal logistic model looking like this:
formula: support_measure ~ effect_mobility + effect_lifequal + effect_health + measure_aq_limit_traffic + effect_smell data: Data_new_reg_sup_meas link threshold nobs logLik AIC niter max.grad cond.H logit flexible 3428 -242.70 521.41 8(0) 1.43e-07 1.9e+02 Coefficients: Estimate Std. Error z value Pr(>|z|) effect_mobilitysehr verbessern 2.1829 0.7348 2.971 0.002972 ** effect_mobilityverbessern 0.9145 0.5250 1.742 0.081536 . effect_mobilityverschlechtern -1.4085 0.3356 -4.197 2.70e-05 *** effect_mobilitysehr verschlechtern -2.7106 0.3517 -7.707 1.29e-14 *** effect_lifequalsehr verbessern 2.9934 0.9742 3.072 0.002123 ** effect_lifequalverbessern 0.9899 0.4044 2.448 0.014368 * effect_lifequalverschlechtern -0.6849 0.4105 -1.668 0.095266 . effect_lifequalsehr verschlechtern -2.0562 0.5684 -3.617 0.000298 *** effect_healthsehr verbessern 0.5033 0.9709 0.518 0.604206 effect_healthverbessern 1.9633 0.4367 4.496 6.93e-06 *** effect_healthverschlechtern 0.3893 0.4400 0.885 0.376355 effect_healthsehr verschlechtern 1.9282 0.5402 3.569 0.000358 *** measure_aq_limit_trafficJa 2.6409 0.2725 9.692 < 2e-16 *** effect_smellsehr verbessern 0.5936 0.9377 0.633 0.526711 effect_smellverbessern 0.5261 0.4588 1.147 0.251501 effect_smellverschlechtern -0.8122 0.3688 -2.202 0.027636 * effect_smellsehr verschlechtern -0.9458 0.4066 -2.326 0.019997 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Threshold coefficients: Estimate Std. Error z value Nein|Ja 2.3023 0.3883 5.93
All of these results make sense, given the distribution of responses in the contingency tables for each predictor vs the response, except for that of
effect_health. When you look at the contingency table of this variable with the response, you see that the estimates don't seem to make any sense.
table(Data_new_reg_sup_meas$support_measure, Data_new_reg_sup_meas$effect_health) kein wirkung sehr verschlechtern verschlechtern verbessern sehr verbessern Nein 931 1215 868 50 7 Ja 56 22 17 162 100
The model estimates indicate that a response of "sehr verschlechtern" (German for 'greatly worsen') is significantly associated with higher odds of the response "Ja" ('yes') for the dependent variable
support_measure. This simply cannot be correct, given the distribution of responses in the contingency table.
When I run the model with just one predictor,
effect_health, then the regression parameters change significantly and seem to accurately describe the underlying relationship:
formula: support_measure ~ effect_health data: Data_new_reg_sup_meas link threshold nobs logLik AIC niter max.grad cond.H logit flexible 3428 -551.20 1112.39 7(0) 7.88e-11 2.9e+01 Coefficients: Estimate Std. Error z value Pr(>|z|) effect_healthsehr verbessern 5.4702 0.4145 13.198 < 2e-16 *** effect_healthverbessern 3.9865 0.2124 18.771 < 2e-16 *** effect_healthverschlechtern -1.1221 0.2809 -3.995 6.48e-05 *** effect_healthsehr verschlechtern -1.2005 0.2554 -4.701 2.58e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Threshold coefficients: Estimate Std. Error z value Nein|Ja 2.8109 0.1376 20.43
My question here is, how is it possible that the regression parameters change so drastically between a model containing all the predictors and one containing one individual predictor?
From a publication by R.H.B. Christensen (cited below), the author of the
ordinal package, I found a tentative explanation:
"Observe that the regression parameter in cumulative link models, cf. (16) are signal-to- noise ratios. This means that adding a covariate to a cumulative link model that reduces the residual noise in the corresponding latent model will increase the signal-to-noise ratios. Thus adding a covariate will (often) increase the coefficients of the other covariates in the cumulative link model. This is different from linear models, where (in orthogonal designs) adding a covariate does not alter the value of the other coefficients"
Christensen R H B 2015 Analysis of ordinal data with cumulative link models—estimation with the ordinal package CRAN R-Project 1–31 Online: https://cran.r-project.org/web/packages/ordinal/vignettes/clm_intro.pdf
While this seems to explain why regression parameters will change when more variables are added to a model, I cannot seem to apply this idea to the major changes I am seeing between these two models for the variable
I have already consulted Alan Agresti's excellent textbooks "Categorical Data Analysis" and "Analysis of Ordinal Categorical Data", as well as an endless set of online resources and have not been able to answer my question. Any responses or pointers of how to answer this question would be greatly appreciated. This is also the first time I'm asking such a question here in Stack Exchange, though I have used it often to help with other issues in the past.
Furthermore, since I am only looking to get the Odds Ratios of the relationships between the predictors and the response variable, does it make more sense to make a model using all of them together and extract the estimates (such as the first example I gave here), or to compare them individually and then extract the estimates (as with the second example I gave)?