#===============================================
# Illustration of cumulative logit model
# using data on mental health from section 6.3.3
#===============================================

# Read mental health data from the web:
data="http://www.stat.ufl.edu/~aa/glm/data/Mental.dat"
mental = read.table(data, header = T)
head(mental)

# We use library VGAM:
library(VGAM)

# We fit a cumulative logit model with main effects of "ses" and "life":
fit.imp=vglm(impair~life+ses,family=cumulative(parallel=T), data=mental)
summary(fit.imp)

# We then consider a model with interaction between "ses" and "life":
fit.int=vglm(impair~life+ses+life:ses,family=cumulative(parallel=T),data=mental)
summary(fit.int)

# We test if there is a significant interaction effect:
G2=deviance(fit.imp)-deviance(fit.int)
df.diff=df.residual(fit.imp)-df.residual(fit.int)
1-pchisq(G2,df.diff)
# The effect of interaction is not significant


# We consider a model where the effects of the covariates may 
# differ between the cumulative logits:
fit.nopar=vglm(impair~life+ses,family=cumulative,data=mental)
summary(fit.nopar)

# We test if this is a significantly better model:
G2=deviance(fit.imp)-deviance(fit.nopar)
df.diff=df.residual(fit.imp)-df.residual(fit.nopar)
1-pchisq(G2,df.diff)