Ordered Logistic Family for Ordinal Regression

Description

Defines the ordered logistic (cumulative logit) family for modeling ordinal response variables within the gamlss2 framework. The location and ordered cutpoints can be modeled as functions of covariates.

Usage

OL(k)

ologit(k)

Arguments

k	A single integer specifying the number of response categories. Must be k >= 2.

Details

This family implements a cumulative logit model for an ordinal response with k categories. The response must be coded as the integers 1, …, k, or as a factor with exactly these levels.

The linear predictor models a latent location parameter. Let \(\gamma_j\) denote the cutpoint between categories \(j\) and \(j + 1\). The cutpoints are parameterized as

\(\gamma_1 = \theta_1\)

and, recursively,

j = {j-1} + (_j), j = 2, , k - 1.

Equivalently,

\(\gamma_j = \gamma_1 + \sum_{\ell=2}^j \exp(\delta_\ell),\)

which ensures strictly ordered cutpoints. In the family object, \(\theta_1\) is named theta1 and \(\delta_j\) is named delta\(j\). For k = 2, only location and theta1 are required; there are no delta parameters.

The cumulative probabilities are

(Y > j) = ^{-1}( - _j), j = 1, , k - 1.

The OL() family supports modeling the location, first cutpoint, and cutpoint differences as functions of covariates using additive predictors in gamlss2 via the “|” formula interface.

The returned family object provides functions for evaluating probability masses and log-likelihoods, cumulative probabilities and quantiles, category and transition probabilities, and the mean and variance of the numeric category labels.

ologit() is retained as a compatibility alias for OL().

Value

A “gamlss2.family” object to be used with gamlss2. Its distribution parameters are location, theta1, and, when k > 2, one cutpoint-difference parameter delta\(j\) for each \(j = 2, \ldots, k - 1\).

See Also

gamlss2, gamlss2.family, polr

Examples

library("gamlss2")

## example using the housing data from the MASS package
library("MASS")

## fit standard cumulative logit model using polr()
m1 <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
summary(m1)

Call:
polr(formula = Sat ~ Infl + Type + Cont, data = housing, weights = Freq)

Coefficients:
                Value Std. Error t value
InflMedium     0.5664    0.10465   5.412
InflHigh       1.2888    0.12716  10.136
TypeApartment -0.5724    0.11924  -4.800
TypeAtrium    -0.3662    0.15517  -2.360
TypeTerrace   -1.0910    0.15149  -7.202
ContHigh       0.3603    0.09554   3.771

Intercepts:
            Value   Std. Error t value
Low|Medium  -0.4961  0.1248    -3.9739
Medium|High  0.6907  0.1255     5.5049

Residual Deviance: 3479.149 
AIC: 3495.149 

## OL() uses integer category labels 1, ..., k
housing$Satint <- as.integer(housing$Sat)

## fit the corresponding model using gamlss2
m2 <- gamlss2(Satint ~ Infl + Type + Cont,
  data = housing, weights = Freq,
  family = OL(k = 3))

GAMLSS-RS iteration  1: Global Deviance = 3479.1765 eps = 0.058036     
GAMLSS-RS iteration  2: Global Deviance = 3479.1493 eps = 0.000007     

summary(m2)

Call:
gamlss2(formula = Satint ~ Infl + Type + Cont, data = housing, 
    family = OL(k = 3), weights = Freq)
---
Family: Ordered Logit (3 categories) 
Link functions: identity, identity, identity
*--------
Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
location.(Intercept)    -0.1969   223.6070  -0.001   0.9993  
location.InflMedium      0.5665     0.5366   1.056   0.2952  
location.InflHigh        1.2889     0.5570   2.314   0.0239 *
location.TypeApartment  -0.5724     0.6343  -0.902   0.3703  
location.TypeAtrium     -0.3662     0.6362  -0.576   0.5669  
location.TypeTerrace    -1.0911     0.6400  -1.705   0.0931 .
location.ContHigh        0.3603     0.4461   0.808   0.4223  
theta1.(Intercept)      -0.6932   223.6070  -0.003   0.9975  
delta2.(Intercept)       0.1714     0.2191   0.782   0.4371  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
*--------
n = 72 df =  9 res.df =  63
Deviance = 3479.1493 Null Dev. Red. = -2099.21%
AIC = 3497.1493 elapsed =  0.03sec

## compare coefficients
coef(m1)

   InflMedium      InflHigh TypeApartment    TypeAtrium   TypeTerrace 
    0.5663937     1.2888191    -0.5723501    -0.3661866    -1.0910149 
     ContHigh 
    0.3602841 

coef(m2)

  location.p.(Intercept)    location.p.InflMedium      location.p.InflHigh 
              -0.1968537                0.5664539                1.2889495 
location.p.TypeApartment    location.p.TypeAtrium   location.p.TypeTerrace 
              -0.5724084               -0.3662225               -1.0911265 
     location.p.ContHigh     theta1.p.(Intercept)     delta2.p.(Intercept) 
               0.3603203               -0.6931558                0.1713817 

## predict class probabilities
pm1 <- predict(m1, type = "p")
pm2 <- family(m2)$probabilities(predict(m2))

print(head(pm1))

        Low    Medium      High
1 0.3784493 0.2876752 0.3338755
2 0.3784493 0.2876752 0.3338755
3 0.3784493 0.2876752 0.3338755
4 0.2568264 0.2742122 0.4689613
5 0.2568264 0.2742122 0.4689613
6 0.2568264 0.2742122 0.4689613

print(head(pm2))

       Pr(Y=1)   Pr(Y=2)   Pr(Y=3)
[1,] 0.3784101 0.2876996 0.3338904
[2,] 0.3784101 0.2876996 0.3338904
[3,] 0.3784101 0.2876996 0.3338904
[4,] 0.2567831 0.2742239 0.4689929
[5,] 0.2567831 0.2742239 0.4689929
[6,] 0.2567831 0.2742239 0.4689929