library(foreign)
library(grid)
library(gridExtra)

cc <- read.table("https://teaching.sociology.ul.ie/so5032/creditcard.dat", col.names=c("income", "card"))

aggregate(cc$income ~ cc$card, FUN=mean)

library(ggplot2)
ggplot(data=cc, aes(x=income, y=card)) + geom_point(color="red", alpha=0.2)
ggsave("credit.png", units="cm", width=20, height=11)

summary(modlin <- lm(data=cc, card ~ income))

cc$plin <- predict(modlin)
ggplot(data=cc, aes(x=income, y=card)) +
    geom_point(color="red", alpha=0.2) +
    geom_point(color="darkorange",aes(x=income,y=plin))
ggsave("creditLPMR.png", units="cm", width=20, height=11)

cc$rlin <- cc$card - cc$plin
p1 <- ggplot(data=cc, aes(x=rlin)) + geom_histogram()
p2 <- ggplot(data=cc, aes(x=income, y=rlin)) + geom_point(color="red", alpha=0.4)
grobj <- arrangeGrob(p1, p2, ncol=2)
ggsave("creditLPMresR.png", grobj, units="cm", width=20, height=11)

ggplot(data.frame(x=c(0, 0.975)), aes(x)) + stat_function(fun=function(x) (x/(1-x)), n=1001) +
  xlab("Probability") + ylab("Odds")
ggsave("p2oddsR.png", units="cm", width=20, height=11)

ggplot(data.frame(x=c(0, 1)), aes(x)) + stat_function(fun=function(x) log(x/(1-x)), n=1001) +
  xlab("Probability") + ylab("Log-Odds")
ggsave("p2logoddsR.png", units="cm", width=20, height=11)

ggplot(data.frame(x=c(0, 1)), aes(x)) + stat_function(fun=function(x) log(x/(1-x)), n=1001) +
  xlab("Probability") + ylab("Log-Odds") + coord_flip()
ggsave("logodds2pR.png", units="cm", width=20, height=11)

summary(modlogit <- glm(data=cc, card ~ income, family="binomial"))

cc$plog <- predict(modlogit, type="response")
cc$plograw <- predict(modlogit)
cc$plog2 <- exp(cc$plograw)/(1+exp(cc$plograw))
ggplot(data=cc, aes(x=income, y=card)) +
    geom_point(color="red", alpha=0.2) +
    geom_point(color="darkorange", aes(x=income,y=plin)) +
    geom_point(color="blue", aes(x=income,y=plog))

ggsave("creditlogitR.png", units="cm", width=20, height=11)

ggplot(data=cc, aes(x=income, y=card)) +
  geom_point(color="red", alpha=0.2) +
  geom_point(color="darkorange", aes(x=income,y=plin)) +
  geom_point(color="blue", aes(x=income,y=plog)) +
  geom_function(data=cc, fun=function (x) 1/(1+ exp(-(modlogit$coefficients["(Intercept)"] + modlogit$coefficients["income"]*x))))

ggsave("sigmoidCCR.png", units="cm", width=20, height=11)