Proportional Odds Model for Ordinal Regression

Fits the proportional odds model to a (preferably ordered) factor response.

Usage

propodds(reverse = TRUE, whitespace = FALSE, ynames = FALSE,
   Thresh = NULL, Trev = reverse, Tref = if (Trev) "M" else 1,
   Intercept = NULL)

Arguments

reverse, whitespace

Logical. Fed into arguments of the same name in cumulative.

ynames

See multinomial for information.

Thresh, Trev, Tref, Intercept

Fed into arguments of the same name in cumulative.

Details

The proportional odds model is a special case from the class of cumulative link models. It involves a logit link applied to cumulative probabilities and a strong parallelism assumption. A parallelism assumption means there is less chance of numerical problems because the fitted probabilities will remain between 0 and 1; however the parallelism assumption ought to be checked, e.g., via a likelihood ratio test. This VGAM family function is merely a shortcut for cumulative(reverse = reverse, link = "logit", parallel = TRUE). Please see cumulative for more details on this model.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm, and vgam.

References

See cumulative.

Author

Thomas W. Yee

Warning

No check is made to verify that the response is ordinal if the response is a matrix; see ordered.

Examples

# Fit the proportional odds model, McCullagh and Nelder (1989,p.179)
pneumo <- transform(pneumo, let = log(exposure.time))
(fit <- vglm(cbind(normal, mild, severe) ~ let, propodds, pneumo))
#> 
#> Call:
#> vglm(formula = cbind(normal, mild, severe) ~ let, family = propodds, 
#>     data = pneumo)
#> 
#> 
#> Coefficients:
#> (Intercept):1 (Intercept):2           let 
#>     -9.676093    -10.581725      2.596807 
#> 
#> Degrees of Freedom: 16 Total; 13 Residual
#> Residual deviance: 5.026826 
#> Log-likelihood: -25.09026 
depvar(fit)  # Sample proportions
#>      normal       mild     severe
#> 1 1.0000000 0.00000000 0.00000000
#> 2 0.9444444 0.03703704 0.01851852
#> 3 0.7906977 0.13953488 0.06976744
#> 4 0.7291667 0.10416667 0.16666667
#> 5 0.6274510 0.19607843 0.17647059
#> 6 0.6052632 0.18421053 0.21052632
#> 7 0.4285714 0.21428571 0.35714286
#> 8 0.3636364 0.18181818 0.45454545
weights(fit, type = "prior")  # Number of observations
#>   [,1]
#> 1   98
#> 2   54
#> 3   43
#> 4   48
#> 5   51
#> 6   38
#> 7   28
#> 8   11
coef(fit, matrix = TRUE)
#>             logitlink(P[Y>=2]) logitlink(P[Y>=3])
#> (Intercept)          -9.676093         -10.581725
#> let                   2.596807           2.596807
constraints(fit)  # Constraint matrices
#> $`(Intercept)`
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> 
#> $let
#>      [,1]
#> [1,]    1
#> [2,]    1
#> 
summary(fit)
#> 
#> Call:
#> vglm(formula = cbind(normal, mild, severe) ~ let, family = propodds, 
#>     data = pneumo)
#> 
#> Coefficients: 
#>               Estimate Std. Error z value Pr(>|z|)    
#> (Intercept):1  -9.6761     1.3241  -7.308 2.72e-13 ***
#> (Intercept):2 -10.5817     1.3454  -7.865 3.69e-15 ***
#> let             2.5968     0.3811   6.814 9.50e-12 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3])
#> 
#> Residual deviance: 5.0268 on 13 degrees of freedom
#> 
#> Log-likelihood: -25.0903 on 13 degrees of freedom
#> 
#> Number of Fisher scoring iterations: 4 
#> 
#> 
#> Exponentiated coefficients:
#>      let 
#> 13.42081 

# Check that the model is linear in let ----------------------
fit2 <- vgam(cbind(normal, mild, severe) ~ s(let, df = 2), propodds,
             pneumo)
if (FALSE)  plot(fit2, se = TRUE, lcol = 2, scol = 2)  # \dontrun{}

# Check the proportional odds assumption with a LRT ----------
(fit3 <- vglm(cbind(normal, mild, severe) ~ let,
              cumulative(parallel = FALSE, reverse = TRUE), pneumo))
#> 
#> Call:
#> vglm(formula = cbind(normal, mild, severe) ~ let, family = cumulative(parallel = FALSE, 
#>     reverse = TRUE), data = pneumo)
#> 
#> 
#> Coefficients:
#> (Intercept):1 (Intercept):2         let:1         let:2 
#>     -9.593308    -11.104791      2.571300      2.743550 
#> 
#> Degrees of Freedom: 16 Total; 12 Residual
#> Residual deviance: 4.884404 
#> Log-likelihood: -25.01905 
pchisq(deviance(fit) - deviance(fit3),
       df = df.residual(fit) - df.residual(fit3), lower.tail = FALSE)
#> [1] 0.7058849
lrtest(fit3, fit)  # Easier
#> Likelihood ratio test
#> 
#> Model 1: cbind(normal, mild, severe) ~ let
#> Model 2: cbind(normal, mild, severe) ~ let
#>   #Df  LogLik Df  Chisq Pr(>Chisq)
#> 1  12 -25.019                     
#> 2  13 -25.090  1 0.1424     0.7059