% \VignetteIndexEntry{Knee Injuries - Proportional Odds Models}
% \VignetteDepends{rms, MASS, VGAM}
%\VignetteEngine{knitr::knitr}
%\VignetteEncoding{UTF-8}
\documentclass[a4paper]{article}

\title{Knee Injuries - Proportional Odds Models}

\begin{document}

\maketitle

<<echo=FALSE,eval=FALSE>>=
options(width=80)
@


At the beginning the knee dataset is loaded:

<<results='hide',eval=FALSE>>=
rm(list=ls(all=TRUE))
library(catdata)
data(knee)
attach(knee)
@

First of all a simple $\chi^{2}$ -- test of independence between the therapy (Th) and the pain level (R4) is applied.
<<eval=FALSE>>=
suppressWarnings(chisq.test(knee$Th,knee$R4))
@

In the following the variable Age is centered around $30$ years and the quadratic variable Age2 is created.
<<eval=FALSE>>=
Age <- Age - 30
Age2 <- Age^2
@

The response pain (R4) has to be an ordered factor, the covariates therapy (Th) and sex (Sex) need to be factors.
<<eval=FALSE>>=
R4 <- as.ordered(R4)
Th <- as.factor(Th)
Sex <- as.factor(Sex)
@

A proportional odds model can be fitted by the function "polr" from the "MASS" -- library. Attention has to be paid to the algebraic signs of the coefficients.
These are inverse to the usual interpretation in porportional odds models.
<<eval=FALSE>>=
library(MASS)
@

The first model only uses therapy as covariate, to achieve a proportional odds model the option "method" needs to use the logistic link function.

<<results='hide',eval=FALSE>>=
polr1 <- polr(R4 ~ Th, method="logistic")
@

<<eval=FALSE>>=
summary(polr1)
@

The corresponding odds--ratio can be recieved by the following command (consider the inverse sign!):
<<eval=FALSE>>=
exp(-coef(polr1))
@



Now a model with the covariates therapy, sex and age is fitted.
<<results='hide',eval=FALSE>>=
polr2 <- polr(R4 ~ Th + Sex + Age, method="logistic")
@

<<eval=FALSE>>=
summary(polr2)
@

Odds--ratios for the second model:
<<eval=FALSE>>=
exp(-coef(polr2))
@

To get the Wald--statistic, the standard errors have to be extracted from the summary.
Afterwards the Wald--statistic and the corresponding p--values are easily recieved.

<<eval=FALSE>>=
se <- summary(polr2)[1][[1]][1:3,2]
(wald2 <- -coef(polr2)/se)
@

P--values for the second model:
<<eval=FALSE>>=
1-pchisq(wald2^2, df=1)
@


Finally the quadratic age--effect is added to the previous model.
<<results='hide',eval=FALSE>>=
polr3 <- update(polr2, ~. + Age2)
@

<<eval=FALSE>>=
summary(polr3)
@

Odds--ratios for the final model:
<<eval=FALSE>>=
exp(-coef(polr3))
@

Wald--statistic for the final model:

<<eval=FALSE>>=
se <- summary(polr3)[1][[1]][1:4,2]
(wald3 <- -coef(polr3)/se)
@

P--values for the final model:
<<eval=FALSE>>=
1-pchisq(wald3^2, df=1)
@

As the proportional odds--model is the most popular model for ordinal data, there are several
different ways to fit such models. Now the final model is additionally fitted with function "vglm" 
from the "VGAM"--library and with function "lrm" from the "rms"--library.

\bigskip

Model fitted with "vglm":

<<eval=FALSE>>=
library(VGAM)
@

<<eval=FALSE>>=
m.vglm <- vglm(R4 ~ Th + Sex + Age + Age2, family = cumulative (link="logit", 
parallel=TRUE))
summary(m.vglm)
@

The resulting coefficients are very similar to the coefficients in the model above fitted
with function "polr", but they have inverse signs. Therefore they can be interpreted in the usual way.

\bigskip

Model fitted with "lrm":
<<eval=FALSE>>=
library(rms)
@

<<eval=FALSE>>=
m.lrm <- lrm(R4 ~ Th + Sex + Age + Age2)
m.lrm
@

Again no big differences are found concerning the coefficients. Here the signs are the same as with the function "polr".

<<echo=FALSE,eval=FALSE>>=
detach(package:VGAM)
@
\end{document}