% \VignetteIndexEntry{Retinopathy - Sequential Logit Models}
% \VignetteDepends{VGAM}
%\VignetteEngine{knitr::knitr}
%\VignetteEncoding{UTF-8}

\documentclass[a4paper]{article}

\title{Retinopathy - Sequential Logit Models}

\begin{document}

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


\maketitle

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

<<results='hide'>>=
library(catdata)
data(retinopathy)
@

For sequential models again the "vglm"--function from the "VGAM"--library is needed, but now family option "sratio" is required.
<<>>=
library(VGAM)
@

Now several sequential logit models are fitted and compared by their corresponding deviances. The first model
is the sequential logit model with all category--specific effects, so the option "parallel=FALSE" is used.

<<>>=
seqm1 <- vglm(RET ~ SM + DIAB + GH + BP, family = sratio (link="logitlink", 
parallel=FALSE), data = retinopathy)
deviance(seqm1)
@

No category--specific effect for DIAB:

<<>>=
seqm2 <- vglm(RET ~ SM + DIAB + GH + BP, family = sratio (link="logitlink", 
parallel=FALSE ~ SM + GH + BP), data = retinopathy)
deviance(seqm2)
@

Testing the removed effect:
<<>>=
1-pchisq(deviance(seqm2)-deviance(seqm1), df=1)
@

No category--specific effect for GH:

<<>>=
seqm3 <- vglm(RET ~ SM + DIAB + GH + BP, family = sratio (link="logitlink", 
parallel=FALSE ~ SM + BP), data = retinopathy)
deviance(seqm3)
@

Testing the removed effect:
<<>>=
1-pchisq(deviance(seqm3)-deviance(seqm2), df=1)
@

No category--specific effect for BP:

<<>>=
seqm4 <- vglm(RET ~ SM + DIAB + GH + BP, family = sratio (link="logitlink", 
parallel=FALSE ~ SM), data = retinopathy)
deviance(seqm4)
@

Testing the removed effect:
<<>>=
1-pchisq(deviance(seqm4)-deviance(seqm3), df=1)
@

No category--specific effect for GH (only global effects):

<<>>=
seqm5 <- vglm(RET ~ SM + DIAB + GH + BP, family = sratio (link="logitlink", 
parallel=TRUE), data = retinopathy)
deviance(seqm5)
@

Testing the removed effect:
<<>>=
1-pchisq(deviance(seqm5)-deviance(seqm4), df=1)
@

As the last test is significant, model "seqm4" is analyzed in detail.

<<>>=
summary(seqm4)
@

The summary gives no p--values for the individual covariates, they have to be computed separately. 
For this purpose the t--values are copied from the summary. The quadratic t--values are the wald--statistics which can be used to produce the individual 
p--values.

\bigskip

p--value intercept1:
<<>>=
1 - pchisq(9.5223^2, df=1)
@

p--value intercept2:
<<>>=
1 - pchisq(8.9957^2, df=1)
@

p--value SM1:
<<>>=
1 - pchisq((-1.8646)^2, df=1)
@

p--value SM2:
<<>>=
1 - pchisq(1.5687^2, df=1)
@

p--value DIAB:
<<>>=
1 - pchisq((-10.4303)^2, df=1)
@

p-value GH:
<<>>=
1 - pchisq((-6.3116)^2, df=1)
@

p--value BP:
<<>>=
1 - pchisq((-5.1037)^2, df=1)
@

To receive the corresponding odds--ratios, the following command can be used.

<<>>=
exp(coefficients(seqm4)[3:7])
@


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