3
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ 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 ... $\endgroup$ – Ben Bolker Dec 27 '16 at 22:32
  • $\begingroup$ @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 $\endgroup$ – Algorithmatic Dec 27 '16 at 23:12
  • $\begingroup$ That's not an easy conclusion to draw in this context. Can you say more about the background and goals of this analysis, please? $\endgroup$ – Ben Bolker Dec 28 '16 at 2:01
4
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ How would you interpret the result of anova in this case? $\endgroup$ – Algorithmatic Dec 29 '16 at 20:29

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.