% \VignetteIndexEntry{Retinopathy - Testing Proportional Odds Assumption}
% \VignetteDepends{VGAM}
%\VignetteEngine{knitr::knitr}
%\VignetteEncoding{UTF-8}
    \documentclass[a4paper]{article}

\title{Retinopathy - Testing Proportional Odds Assumption}

\begin{document}

<<include=FALSE>>=
library(knitr)
opts_chunk$set(
concordance=TRUE
)
@


\maketitle

<<echo=FALSE>>=
options(width=60)
rm(list=ls(all=TRUE))
@

<<results='hide'>>=
library(catdata)
data(retinopathy)
retinopathy$SM <- as.factor(retinopathy$SM)
@

For the fitting of the partial proportional odds models the function "vglm" from the "VGAM"--package is used. First a simple
proportional odds model is fitted with "vglm". 

For the "vglm"--function the response (RET) does not necessarily have to be ordered, SM has to be factorized.

<<>>=
library(VGAM)
# retinopathy$RET <- as.ordered(retinopathy$RET)
# retinopathy$SM <- as.factor(retinopathy$SM)
@


The models differ in the option "parallel" for the used family "cumulative". 

<<results='hide'>>=
pom <- vglm(RET ~ SM + DIAB + GH + BP, family = cumulative(parallel=TRUE), 
            data = retinopathy)
@

<<results='hide'>>=
ppom <- vglm(RET ~ SM + DIAB + GH + BP, family = cumulative(parallel=FALSE), 
             data = retinopathy)
@

First the proportional odds assumption is tested. The deviances of the two models can be received by the following command.

<<>>=
deviance(pom)
@

<<>>=
deviance(ppom)
@

The p--value for the proportional odds assuption is computed:

<<>>=
1 - pchisq(deviance(pom) - deviance(ppom), df=4)
@

Coefficients and standard errors of both models are obtainen in the corresponding summaries. 
\bigskip

Summary proportional odds model:
<<>>=
summary(pom)
@

\bigskip
Summary partial proportional odds model:
<<>>=
summary(ppom)
@

\bigskip

Now the proportional odds assumption for all covariates is taken away step by step. 
Afterwards the corresponding proportional odds assumptions are tested. 

\bigskip
 
Global effect for BP:
<<>>=
ppom2 <- vglm (RET ~ SM + DIAB + GH + BP,
family = cumulative (parallel = FALSE ~ SM + DIAB + GH), data = retinopathy)
deviance(ppom2)
@

<<>>=
1-pchisq(deviance(ppom2)-deviance(ppom), df=1)
@


Global effect for GH:
<<>>=
ppom3 <- vglm (RET ~ SM + DIAB + GH + BP, 
family = cumulative (parallel = FALSE ~ SM + DIAB), data = retinopathy)
deviance(ppom3)
@

<<>>=
1-pchisq(deviance(ppom3)-deviance(ppom2), df=1)
@

Global effect for DIAB:
<<>>=
ppom4 <- vglm (RET ~ SM + DIAB + GH + BP, 
family = cumulative (parallel = FALSE ~ SM), data = retinopathy)
deviance(ppom4)
@

<<>>=
1-pchisq(deviance(ppom4)-deviance(ppom3), df=1)
@

Global effect for SM (equivalent to proportional odds model):
<<>>=
1-pchisq(deviance(pom)-deviance(ppom4), df=1)
@

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