# GLM interpretation of parameters of ordinal predictor variables

When using the glm function in the stats package with ordinal variables

data$X1 <- ordered(data$X1)
data$X2 <- ordered(data$X2)
data$X3 <- ordered(data$X3)
data$X4 <- ordered(data$X4)
data$X5 <- ordered(data$X5)
data$X6 <- ordered(data$X6)

glm(Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = data)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  2.19445    0.07169  30.610  < 2e-16 ***
X1.L         0.48235    0.10060   4.795 2.48e-06 ***
X1.Q        -0.42570    0.08366  -5.089 6.11e-07 ***
X1.C        -0.01543    0.07109  -0.217  0.82832
X1^4        -0.06392    0.05936  -1.077  0.28230
X2.L         0.51696    0.26341   1.963  0.05055 .
X2.Q        -0.34153    0.19369  -1.763  0.07879 .
X2.C         0.28979    0.13800   2.100  0.03651 *
X2^4        -0.30389    0.10082  -3.014  0.00278 **
X3.L        -0.01790    0.10756  -0.166  0.86790
X3.Q         0.17086    0.08730   1.957  0.05119 .
X3.C         0.04131    0.07452   0.554  0.57975
X3^4         0.08683    0.06149   1.412  0.15891
X4.L         0.18349    0.27466   0.668  0.50457
X4.Q         0.32422    0.20563   1.577  0.11583
X4.C        -0.39733    0.14880  -2.670  0.00796 **
X4^4         0.29617    0.10123   2.926  0.00368 **
X5.L         0.29977    0.10860   2.760  0.00610 **
X5.Q        -0.35389    0.08825  -4.010 7.53e-05 ***
X5.C         0.03904    0.07258   0.538  0.59102
X5^4        -0.02931    0.05908  -0.496  0.62011
X6.L         0.16479    0.18020   0.915  0.36113
X6.Q        -0.14659    0.15071  -0.973  0.33146
X6.C         0.05724    0.11289   0.507  0.61245
X6^4        -0.09772    0.07570  -1.291  0.19764
---
...


I understand the variables with .L, .Q, and .C and ^4 are, respectively, the coefficients for the ordered factor coded with linear, quadratic, cubic, and quartic contrasts.

Some of the variables are significant (***, **, *) but they're nonlinear (for example: X2). My question is, when reporting the result of which variable is statically significant in predicting Y, and whether these variable have positive weights or not, which coefficient should you use? These predictors (X1,X2,X4, X5) have significant coefficients in their linear, quadratic, cubic etc. forms, but in opposite directions.

Is it valid to say that X2 is statically significant in predicting Y and has a positive correlation?

Can someone point me to something that would help me understand how to report the results of glm analyses using ordered factors?

• probably to test overall effects of X2 via updating the model (model2 <- update(model1, . ~ . - X2)) and doing a likelihood ratio test (anova(model1, model2)). Would be hard to make a statement about overall positive correlations in a model like this ... Dec 27 '16 at 22:32
• @BenBolker My understanding is that doing it this way would test which model (model1 without X2 and model2 with X2) offers significantly better goodness-of-fit, but doesn't indicate if X2 has a positive correlation with Y Dec 27 '16 at 23:12
• That's not an easy conclusion to draw in this context. Can you say more about the background and goals of this analysis, please? Dec 28 '16 at 2:01

I'm using a slightly cooked example here.

data("kyphosis",package="HH")
kk <- subset(kyphosis,Number<=5 & Start>=13 & Start<=17)
kk <- transform(kk, Number=ordered(Number), Start=ordered(Start),
kyph=as.numeric(Kyphosis)-1)


Unfortunately since I've cut the data set down so much there are only 2 instances of kyphosis presnt (out of 43 cases), so I have to use the brglm package to overcome complete separation

library(brglm)
m1 <- brglm(Kyphosis~Number+Start,family=binomial,data=kk)

## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.038967   0.581733  -3.505 0.000457 ***
## Number.L     0.221716   1.321268   0.168 0.866736
## Number.Q     0.074500   1.286798   0.058 0.953832
## Number.C     0.359261   1.227672   0.293 0.769801
## Start.L     -0.813621   1.627716  -0.500 0.617178
## Start.Q      0.295470   1.584274   0.187 0.852051
## Start.C      1.231051   1.259250   0.978 0.328269
## Start^4      0.007603   1.430817   0.005 0.995760


The test:

 m2 <- update(m1, . ~ . - Number)
anova(m1,m2,test="Chisq")


This tests whether Number has a significant effect overall on the incidence of kyphosis. Determining whether Number significantly increases the incidence of kyphosis is tricky, because the higher-level/nonlinear contrasts (quadratic/Q, cubic/C, quartic/4, etc.) mean that the effect of increasing Number by 1 unit can depend on where you start. For example, if the quadratic term is large and (let's say, without much loss of generality) positive, then the effect of increasing Number will generally be to decrease the incidence when Number is small and to increase it when Number is large. I suppose that if the linear term is positive and significant and all the other terms are small, then you could say there is a positive effect of Number on kyphosis, but otherwise it will be difficult to say.

• How would you interpret the result of anova in this case? Dec 29 '16 at 20:29