Ordinal regression

TODO

• link to regressionLogistic, regressionMultinom

Install required packages

MASS, ordinal, rms, VGAM, DescTools

wants <- c("MASS", "ordinal", "rms", "VGAM", "DescTools")
has   <- wants %in% rownames(installed.packages())
if(any(!has)) install.packages(wants[!has])

Ordinal regression (proportional odds model)

Simulate data

Dependent variable \(Y_{\text{ord}}\) with \(k=4\) groups, \(p=2\) predictor variables

set.seed(123)
N     <- 100
X1    <- rnorm(N, 175, 7)
X2    <- rnorm(N,  30, 8)
Ycont <- 0.5*X1 - 0.3*X2 + 10 + rnorm(N, 0, 6)
Yord  <- cut(Ycont, breaks=quantile(Ycont), include.lowest=TRUE,
             labels=c("--", "-", "+", "++"), ordered=TRUE)
dfOrd <- data.frame(X1, X2, Yord)

Using vglm() from package VGAM

Model using cumulative logits: \(\text{logit}(p(Y \geq g)) = \ln \frac{P(Y \geq g)}{1 - P(Y \geq g)} = \beta_{0_{g}} + \beta_{1} X_{1} + \dots + \beta_{p} X_{p} \quad(g = 2, \ldots, k)\)

library(VGAM)
(vglmFit <- vglm(Yord ~ X1 + X2, family=propodds, data=dfOrd, model=TRUE))

Call:
vglm(formula = Yord ~ X1 + X2, family = propodds, data = dfOrd, 
    model = TRUE)


Coefficients:
(Intercept):1 (Intercept):2 (Intercept):3            X1            X2 
 -15.61123204  -17.00112492  -18.28506734    0.11197395   -0.09517965 

Degrees of Freedom: 300 Total; 295 Residual
Residual deviance: 249.3579 
Log-likelihood: -124.6789 

Equivalent:

vglm(Yord ~ X1 + X2, family=cumulative(parallel=TRUE, reverse=TRUE), data=dfOrd)
# not shown

Adjacent category logits \(\ln \frac{P(Y=g)}{P(Y=g-1)}\) with proportional odds assumption

vglm(Yord ~ X1 + X2, family=acat(parallel=TRUE), data=dfOrd)
# not shown

Continuation ratio logits \(\ln \frac{P(Y=g)}{P(Y < g)}\) with proportional odds assumption (discrete version of Cox proportional hazards model for survival data)

vglm(Yord ~ X1 + X2, family=sratio(parallel=TRUE), data=dfOrd)
# not shown

Using orm() from package rms

Model \(\text{logit}(p(Y \geq g)) = \beta_{0_{g}} + \beta_{1} X_{1} + \dots + \beta_{p} X_{p} \quad(g = 2, \ldots, k)\)

library(rms)
(ormFit <- orm(Yord ~ X1 + X2, data=dfOrd))

Logistic (Proportional Odds) Ordinal Regression Model
 
 orm(formula = Yord ~ X1 + X2, data = dfOrd)
 
 
 Frequencies of Responses
 
 --  -  + ++ 
 25 25 25 25 
 
                       Model Likelihood          Discrimination          Rank Discrim.    
                          Ratio Test                 Indexes                Indexes       
 Obs           100    LR chi2      27.90    R2                  0.260    rho     0.477    
 Distinct Y      4    d.f.             2    g                   1.176                     
 Median Y        2    Pr(> chi2) <0.0001    gr                  3.240                     
 max |deriv| 0.003    Score chi2   28.50    |Pr(Y>=median)-0.5| 0.274                     
                      Pr(> chi2) <0.0001                                                  
 
       Coef     S.E.   Wald Z Pr(>|Z|)
 y>=-  -15.6110 5.5109 -2.83  0.0046  
 y>=+  -17.0008 5.5508 -3.06  0.0022  
 y>=++ -18.2848 5.5863 -3.27  0.0011  
 X1      0.1120 0.0314  3.56  0.0004  
 X2     -0.0952 0.0272 -3.50  0.0005  
 

Using polr() from package MASS

Model \(\text{logit}(p(Y \leq g)) = \beta_{0_{g}} - (\beta_{1} X_{1} + \dots + \beta_{p} X_{p}) \quad(g = 1, \ldots, k-1)\)

library(MASS)
(polrFit <- polr(Yord ~ X1 + X2, method="logistic", data=dfOrd))
# not shown

Profile likelihood based confidence intervals (need to use MASS:::confint.polr() instead of confint() since other packages are loaded, and method is masked).

exp(MASS:::confint.polr(polrFit))

       2.5 %    97.5 %
X1 1.0530865 1.1919021
X2 0.8602671 0.9574481

Using clm() from package ordinal

Model \(\text{logit}(p(Y \leq g)) = \beta_{0_{g}} - (\beta_{1} X_{1} + \dots + \beta_{p} X_{p}) \quad(g = 1, \ldots, k-1)\)

library(ordinal)
(clmFit <- clm(Yord ~ X1 + X2, link="logit", data=dfOrd))
# not shown

Predicted category membership

Predicted category probabilities

PhatCateg <- VGAM::predict(vglmFit, type="response")
head(PhatCateg)

          --         -         +        ++
1 0.22610471 0.3136747 0.2692008 0.1910199
2 0.32021125 0.3338845 0.2181580 0.1277463
3 0.07320949 0.1675519 0.2930451 0.4661935
4 0.19019915 0.2950991 0.2876648 0.2270369
5 0.12403581 0.2383874 0.3099813 0.3275955
6 0.07534083 0.1711326 0.2950389 0.4584877

predict(ormFit, type="fitted.ind")
predict(clmFit, subset(dfOrd, select=c("X1", "X2"), type="prob"))$fit
predict(polrFit, type="probs")
# not shown

Predicted categories

categHat <- levels(dfOrd$Yord)[max.col(PhatCateg)]
head(categHat)

[1] "-"  "-"  "++" "-"  "++" "++"

predict(clmFit, type="class")
predict(polrFit, type="class")
# not shown

Apply regression model to new data

Simulate new data

Nnew  <- 3
dfNew <- data.frame(X1=rnorm(Nnew, 175, 7),
                    X2=rnorm(Nnew,  30, 8))

Predicted class probabilities

VGAM::predict(vglmFit, dfNew, type="response")

         --          -          +          ++
1 0.8809481 0.08648444 0.02333061 0.009236851
2 0.4245288 0.32303923 0.16691358 0.085518387
3 0.1791804 0.28786310 0.29281886 0.240137663

predict(ormFit,  dfNew, type="fitted.ind")
predict(polrFit, dfNew, type="probs")
predict(clmFit,  subset(dfNew, select=c("X1", "X2"), type="prob"))$fit
# not shown

Assess model fit

Classification table

facHat <- factor(categHat, levels=levels(dfOrd$Yord))
cTab   <- xtabs(~ Yord + facHat, data=dfOrd)
addmargins(cTab)

     facHat
Yord   --   -   +  ++ Sum
  --   17   4   3   1  25
  -     5  11   2   7  25
  +     1  10   4  10  25
  ++    3   9   2  11  25
  Sum  26  34  11  29 100

Correct classification rate

(CCR <- sum(diag(cTab)) / sum(cTab))

[1] 0.43

Deviance, log-likelihood and AIC

VGAM::deviance(vglmFit)

[1] 249.3579

VGAM::logLik(vglmFit)

[1] -124.6789

VGAM::AIC(vglmFit)

[1] 259.3579

McFadden, Cox & Snell and Nagelkerke pseudo \(R^{2}\)

library(DescTools)
PseudoR2(vglmFit, which=c("McFadden", "CoxSnell", "Nagelkerke"))

  McFadden   CoxSnell Nagelkerke 
 0.1006315  0.2434676  0.2596987 

Coefficient tests and overall model test

Individual coefficient tests

Estimated standard deviations, z-values and p-values for parameters based on assumption that z-values are asymptotically \(N(0, 1)\) distributed.

sumOrd   <- summary(vglmFit)
(coefOrd <- coef(sumOrd))

                  Estimate Std. Error   z value     Pr(>|z|)
(Intercept):1 -15.61123204 5.41912617 -2.880766 0.0039671060
(Intercept):2 -17.00112492 5.45613579 -3.115964 0.0018334440
(Intercept):3 -18.28506734 5.49803759 -3.325744 0.0008818278
X1              0.11197395 0.03122493  3.586043 0.0003357330
X2             -0.09517965 0.02694012 -3.533007 0.0004108612

Profile likelihood confidence intervals

VGAM::confint(vglmFit, method="profile")

                     2.5 %      97.5 %
(Intercept):1 -26.71111227 -4.99528461
(Intercept):2 -28.19013553 -6.32008987
(Intercept):3 -29.55595252 -7.54783873
X1              0.05172567  0.17554897
X2             -0.15051017 -0.04348481

Tests for other models.

summary(polrFit)
summary(clmFit)
# not shown

Model comparisons - likelihood-ratio tests

Likelihood-ratio-test for predictor X2

vglmR <- vglm(Yord ~ X1, family=propodds, data=dfOrd)
anova(vglmFit, vglmR, type="I")

Analysis of Deviance Table

Model 1: Yord ~ X1 + X2
Model 2: Yord ~ X1
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1       295     249.36                          
2       296     262.84 -1  -13.482 0.0002408 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Likelihood-ratio-test for the full model against the 0-model without predictors (just intercept)

vglm0 <- vglm(Yord ~ 1, family=propodds, data=dfOrd)
anova(vglmFit, vglm0, type="I")

Analysis of Deviance Table

Model 1: Yord ~ X1 + X2
Model 2: Yord ~ 1
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1       295     249.36                          
2       297     277.26 -2  -27.901 8.737e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Test assumption of proportional odds (=parallel logits)

Using vglm() from package VGAM

vglmP  <- vglm(Yord ~ X1 + X2, family=cumulative(parallel=TRUE,  reverse=TRUE),
              data=dfOrd)

# vglmNP <- vglm(Yord ~ X1 + X2, family=cumulative(parallel=FALSE, reverse=TRUE),
#                data=dfOrd)
# VGAM::lrtest(vglmP, vglmNP)

Using clm() from package ordinal

clmP  <- clm(Yord ~ X1 + X2, link="logit", data=dfOrd)

## model with non-proportional odds for X2:
clmNP <- clm(Yord ~ X1, nominal=~X2, data=dfOrd)
anova(clmP, clmNP)

Likelihood ratio tests of cumulative link models:
 
      formula:       nominal: link: threshold:
clmP  Yord ~ X1 + X2 ~1       logit flexible  
clmNP Yord ~ X1      ~X2      logit flexible  

      no.par    AIC  logLik LR.stat df Pr(>Chisq)
clmP       5 259.36 -124.68                      
clmNP      7 259.96 -122.98   3.398  2     0.1829

Detach (automatically) loaded packages (if possible)

try(detach(package:VGAM))
try(detach(package:splines))
try(detach(package:stats4))
try(detach(package:rms))
try(detach(package:Hmisc))
try(detach(package:lattice))
try(detach(package:survival))
try(detach(package:Formula))
try(detach(package:ggplot2))
try(detach(package:MASS))
try(detach(package:ordinal))
try(detach(package:DescTools))

Get the article source from GitHub

R markdown - markdown - R code - all posts