myfile =
happy <- read.table("http://sites.williams.edu/bklingen/files/2013/12/happiness.dat", header=TRUE)
head(happy)
## VGAM Package
require(VGAM)
fit <- vglm(cbind(very,pretty,not) ~ income, data=happy, family=cumulative(parallel=TRUE))
summary(fit)
fit1 <- vglm(cbind(very,pretty,not) ~ factor(income), data=happy, family=cumulative(parallel=TRUE))
summary(fit1)
LR <- -2*(logLik(fit)-logLik(fit1))
1-pchisq(LR,df=1)

### fitted logits
fit.logit <- predict(fit)
name <- colnames(fit.logit)
attach(happy)
plot(fit.logit[,1]~income, type="b", col="red", ylim=c(-2,4.5), ylab="Fitted Cumulative Logit", xlab="Family Income", main="Proportional odds model \n with income as quantitative")
lines(fit.logit[,2]~income, type="b", col="blue")
## add sample logits:
n <- rowSums(happy[,2:4]) #total sample size
sample.cumprob1 <- happy[,2]/n
sample.cumprob2 <- rowSums(happy[,2:3])/n
sample.logit1 <- logit(sample.cumprob1)
sample.logit2 <- logit(sample.cumprob2)
points(sample.logit1~income, pch="+", col="red")
points(sample.logit2~income, pch="+", col="blue")
legend("top", legend=name, lty=c(1,1), col=c("red","blue"), ncol=2, bty=n)

fit1.logit <- predict(fit1)
name <- colnames(fit1.logit)
plot(fit1.logit[,1]~income, type="b", col="red", ylim=c(-2,4.5), ylab="Fitted Cumulative Logit", xlab="Family Income", main="Proportional odds model \n with income as qualitative")
lines(fit1.logit[,2]~income, type="b", col="blue")
points(sample.logit1~income, pch="+", col="red")
points(sample.logit2~income, pch="+", col="blue")
legend("top", legend=name, lty=c(1,1), col=c("red","blue"), ncol=2, bty=n)

### fitted cumulative Probs
fit.cumProb <- predict(fit, untransform=TRUE)
name <- colnames(fit.cumProb)
plot(fit.cumProb[,1]~income, type="b", col="red", lwd=2, ylim=c(0,1), ylab="Fitted Cumulative Probability", xlab="Family Income", main="Proportional odds model")
lines(fit.cumProb[,2]~income, type="b", col="blue", lwd=2)
legend("bottomright", legend=name,lty=c(1,1),col=c("red","blue"), lwd=c(2,2), ncol=1, bty=n)

### fitted category Probs
fit.prob <- predict(fit, type="response")
name <- colnames(fit.prob)
plot(fit.prob[,1]~income, type="b", col="red", lwd=2, ylim=c(0,1), ylab="Fitted Probability", xlab="Family Income", main="Proportional odds model")
lines(fit.prob[,2]~income, type="b", col="blue", lwd=2)
lines(fit.prob[,3]~income, type="b", col="green", lwd=2)
legend("topright", legend=name,lty=c(1,1,1),col=c("red","blue","green"), lwd=c(2,2,2), ncol=1)



#### Non-proportional Odds Model #######
fit2 <- vglm(cbind(very,pretty,not) ~ income, data=happy, family=cumulative(parallel=FALSE))
summary(fit2)
### fitted logits
fit2.logit <- predict(fit2)
name <- colnames(fit2.logit)
plot(fit2.logit[,1]~income, type="b", col="red", ylim=c(-2,4.5), ylab="Fitted Cumulative Logit", xlab="Family Income", main="Model without proportional odds \n with income as quantitative")
lines(fit2.logit[,2]~income, type="b", col="blue")
## add sample logits:
points(sample.logit1~income, pch="+", col="red")
points(sample.logit2~income, pch="+", col="blue")
legend("top", legend=name,lty=c(1,1),col=c("red","blue"), ncol=2, bty=n)


#### Adjacent category logit model   #######
# by default, VGLM fits log(P[Y=j+1]/P[Y=j]), reverse=TRUE reverses this to log(P[Y=j]/P[Y=j+1])
fit3 <- vglm(cbind(very,pretty,not) ~ income, ,data=happy, family=acat(reverse=TRUE,parallel=TRUE))
summary(fit3)
fitted(fit3)
fit4 <- vglm(cbind(very,pretty,not) ~ income, data=happy, family = acat(reverse=TRUE))
summary(fit4)

#Goodness of Fit for proportional odds
obs <- happy[2:4]
n <- rowSums(obs)
exp <- apply(fitted(fit2),2,function(col) col*n)
X2 <- sum((obs-exp)^2/exp)
G2 <- 2*sum(obs*log(obs/exp))

# Continuation Ratio logit model
# by default, VGLM fits log(P[Y=j+1]/P[Y=j]), reverse=TRUE reverses this to log(P[Y=j]/P[Y=j+1])
#fit5 <- vglm(cbind(very,pretty,not) ~ income, data=happy, family=cratio(parallel=TRUE))
# or: use sratio() option for family:
fit5 <- vglm(cbind(very,pretty,not) ~ income, data=happy, family=sratio(parallel=TRUE))
summary(fit5)

#### Cumulative Probit Model   ######
fit6 <- vglm(cbind(very,pretty,not) ~ income, data=happy, family=cumulative(link=probit, parallel=TRUE))
summary(fit6)
fit.cumProb.clogit <- predict(fit, untransform=TRUE)
fit.cumProb.cprobit <- predict(fit6, untransform=TRUE)
plot(fit.cumProb.clogit[,1]~income, type="b", lwd=2, col="brown", ylim=c(0,1), ylab="Fitted Cumulative Probability", xlab="Family Income", main="Comparison of cumulative \n logit and probit model")
lines(fit.cumProb.clogit[,2]~income, type="b", lwd=2, col="brown")
lines(fit.cumProb.cprobit[,1]~income, type="b", lwd=2, col="orange")
lines(fit.cumProb.cprobit[,2]~income, type="b", lwd=2, col="orange")
legend("bottomright", legend=c("cum. logit","cum. probit"),lty=c(1,1),lwd=c(2,2), col=c("brown","orange"), ncol=1, bty=n)