happy <- read.csv("http://www.hofroe.net/stat557/data/happy.csv", strip.white=TRUE)
happy[happy=="NA"] <- NA

summary(happy)

library(productplots)

happy$happy <- factor(happy$happy)
happy$age <- as.numeric(as.character(happy$age))

library(MASS)
help(polr)

happy.age <- polr(happy~poly(age,4)*sex, data=na.omit(happy[,c("happy","age","sex")]))

X <- data.frame(expand.grid(
	age=18:80, sex=c("female","male")
))
X <- data.frame(X,predict(happy.age, newdata=X, type="p"))
library(reshape2)
X.melt <- melt(X, id.vars=c("age","sex"))

library(ggplot2)
qplot(age, value, geom=c("point","line"), colour=variable, shape=sex, group=interaction(sex, variable),data=X.melt)+ylim(c(0,1))+ylab("Estimated Probabilities")