3
$\begingroup$

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 effect_health.

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)?

$\endgroup$
  • 1
    $\begingroup$ You're aware this happens with ordinary least squares regression as well? $\endgroup$ – Glen_b Dec 18 '18 at 14:24
  • $\begingroup$ @Glen_b yes I'm aware that this can happen in OLSR, but what concerns me is the magnitude of the changes here. The estimate for the category "sehr verschlechtern" changes from -1.2005 (in the model as a singular predictor) to 1.9282 (in the model with other predictors), which completely changes the direction and strength of the relationship. $\endgroup$ – Sean Schmitz Dec 18 '18 at 15:05
  • 1
    $\begingroup$ Such changes can be quite large - potentially be of any magnitude. $\endgroup$ – Glen_b Dec 18 '18 at 15:27
  • $\begingroup$ @Glen_b would you have any recommendations for my second question then?: 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)? $\endgroup$ – Sean Schmitz Dec 18 '18 at 15:44
  • 1
    $\begingroup$ First consider whether you may have still other important predictors not in the model (also see here). Without a controlled randomized experiment, you can't ignore this possibility. If you don't have any potential omitted variables then the full model should yield unbiased estimates, but if there are many IVs you may want to have some shrinkage. It may help to read about instrumental variables $\endgroup$ – Glen_b Dec 18 '18 at 22:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.