Ordered Logistic Family for Ordinal Regression

Description

Defines the ordered logistic (cumulative logit) family for modeling ordinal response variables using gamlss2. 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 = 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) 
Link functions: identity, identity, identity
*--------
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.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.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