4

I have a dichotomous outcome (gallstones/no gallstones) and an ordinal predictor variable consisting of four classes (body mass index <25(ref.), 25-30, 30-35, 35-45). I want to perform a test for trends using R.

When running the glm function only the bmi 35-45 is significantly associated with the outcome (gallstones):

glm(formula = gallstones ~ bmi, family = binomial, data = m)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max

-0.9282  -0.4058  -0.3624  -0.3624   2.3477


Coefficients:
               Estimate Std. Error z value Pr(>|z|)

(Intercept)     -2.6902     0.1476 -18.223  < 2e-16 ***


bmi(25,30]   0.2347     0.2131   1.102   0.2706

bmi(30,35]   0.6108     0.3291   1.856   0.0635 .

bmi(35,45]   2.0712     0.4915   4.214 2.51e-05 ***

To examine if there is a dose-response relationship with a linear trend with rising bmi I perform the analysis of variance with the chi-square test:

anova(inc_test_3, test="Chisq")
Analysis of Deviance Table

Model: binomial, link: logit


Response: gallstones


Terms added sequentially (first to last)


      Df Deviance Resid. Df Resid. Dev Pr(&gt;Chi)   



NULL                      1489     810.05

bmi       3   15.892      1486     794.16 0.001193 **

There is a significant difference between the obtained estimates from the ordinal bmi predictor.

Is this the correct way to test for a trend on an ordinal variable?

utobi
  • 11,726
Daniel
  • 41
  • 1
    In your analysis of deviance you haven't tested for trend at all, but simply compared the null model with one having a four-level categorical predictor. Orthogonal polynomial coding is useful for examining trends, & in fact is the default when you tell R the predictor is ordinal (see ?ordered). – Scortchi - Reinstate Monica Apr 13 '15 at 11:59
  • Thank you - I was not aware of the ordered function in R. However, what if we concider the variable to be a factor? I am not sure if it is all right to assume that the jumps in this BMI variable are completely ordinal. I have created the BMI variable with the default cut function in R. – Daniel Apr 13 '15 at 13:31
  • (1) Not sure what you mean by "consider the variable to be a factor". Orthogonal polynomial coding doesn't produce a substantively different model from e.g. reference level coding. (2) It doesn't seem at all likely that the chance of having a gallstone jumps at BMIs of 25, 30, 35, & 45, or that it's constant between those BMIs. See here. If you have BMI then use it as a predictor rather than cutting it, allowing for curvilinearity with polynomials or (probably better) splines. – Scortchi - Reinstate Monica Apr 13 '15 at 13:41
  • Sorry for jumping to conclusions to soon:) What i mean is the scenario where the continuous variable has no association on the outcome but I assume it to have an association based on previous analyses. Therefore, I would like to cut the variable into a multi-class categorical variable. Would that be totally wrong? – Daniel Apr 13 '15 at 14:38
  • Sorry - I don't understand the "therefore". Are you saying that when you treated the predictor as having a linear relationship to the response (on the logit scale) you found no (or only a weak) association, & you therefore binned the predictor to allow for a curvilinear relationship? If so, that's not totally wrong, but not a very good way of modelling it - read the post I linked to. – Scortchi - Reinstate Monica Apr 13 '15 at 14:58
  • @Scortchi The orthogonal polynomial coding worked perfect! You want to post it as an answer so i can rate it? – Daniel Apr 18 '15 at 08:42

1 Answers1

2

For the sake of completeness, and expanding over   @Scortchi - Reinstate Monica's comments below the post, the solution is to force R to recognise that bmi is an ordered factor variable. This can be done by calling the ordered function, thus the correct R code would be:

m$ordered_bmi <- ordered(bmi)
glm(formula = gallstones ~ ordered_bmi, family = binomial, data = m)

Ensure that the order is the one you have in mind by checking the levels of odrederd_bmi. If that doesn't happen, feed in the order of the factors yourself by specifying the levels option, e.g.

m$ordered_bmi <- ordered(bmi, levels = c("<25", "25-30", "30-35", "35-45"))
utobi
  • 11,726