Ordered Logistic Family for Ordinal Regression

Description

Defines the ordered logistic (cumulative logit) family for ordinal responses in gamlss2.

Usage

OL(k)

ologit(k)

Arguments

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

Details

The response must be coded as the integers 1, …, k, or as a factor with exactly these levels.

Let \(\gamma_j\) denote the cutpoint between categories \(j\) and \(j + 1\). The cutpoints are parameterized as

\(\gamma_1 = \theta_1\)

and, recursively,

\(\gamma_j = \gamma_{j-1} + \exp(\delta_j), \quad j = 2, \ldots, k - 1.\)

Equivalently,

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

which ensures strictly ordered cutpoints.

The cumulative probabilities are

\(\Pr(Y > j) = \mathrm{logit}^{-1}(\mathrm{location} - \gamma_j), \quad j = 1, \ldots, k - 1.\)

Each parameter can use an additive predictor through the “|” formula interface. The family also provides category and transition probabilities, and the mean and variance of numeric category labels.

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

Value

A “gamlss2.family” object with parameters location, theta1 (\(\theta_1\)), and, for k > 2, delta\(j\) (\(\delta_j\)) for \(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 = 3483.3717 eps = 0.056900     
GAMLSS-RS iteration  2: Global Deviance = 3479.4958 eps = 0.001112     
GAMLSS-RS iteration  3: Global Deviance = 3479.1765 eps = 0.000091     
GAMLSS-RS iteration  4: Global Deviance = 3479.1514 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)(identity(location), identity(theta1),
          identity(delta2)) 
*--------
Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
location.(Intercept)   -0.19494    0.12260  -1.590 0.116833    
location.InflMedium     0.56714    0.10433   5.436 9.41e-07 ***
location.InflHigh       1.29044    0.12589  10.250 4.73e-15 ***
location.TypeApartment -0.57307    0.11904  -4.814 9.63e-06 ***
location.TypeAtrium    -0.36663    0.15521  -2.362 0.021271 *  
location.TypeTerrace   -1.09240    0.15070  -7.249 7.36e-10 ***
location.ContHigh       0.36073    0.09545   3.779 0.000351 ***
theta1.(Intercept)     -0.69316    0.04613 -15.025  < 2e-16 ***
delta2.(Intercept)      0.17236    0.03626   4.753 1.20e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
*--------
n = 72 df =  9 res.df =  63
Deviance = 3479.1514 Null Dev. Red. = -2099.21%
AIC = 3497.1514 elapsed =  0.01sec

## 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.1949367                0.5671416                1.2904360 
location.p.TypeApartment    location.p.TypeAtrium   location.p.TypeTerrace 
              -0.5730733               -0.3666344               -1.0924038 
     location.p.ContHigh     theta1.p.(Intercept)     delta2.p.(Intercept) 
               0.3607349               -0.6931556                0.1723564 

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

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.3779593 0.2879814 0.3340593
[2,] 0.3779593 0.2879814 0.3340593
[3,] 0.3779593 0.2879814 0.3340593
[4,] 0.2562864 0.2743603 0.4693533
[5,] 0.2562864 0.2743603 0.4693533
[6,] 0.2562864 0.2743603 0.4693533