d=read.table("Crabs.txt", header=T)
d$y=factor(d$y,levels=c(0,1),labels=c("No", "Yes"))
d$color=factor(d$color,levels=c(1,2,3,4),labels=c("medium light", "medium", "medium dark", "dark"))
d$spine=factor(d$spine,levels=c(1,2,3),labels=c("both good", "one broker", "both broken"))
str(d)
table(d$y)
hist(d$sat)

##Question 2
fit.null=glm(y~1, data=d, family=binomial)

summary(glm(y~weight, data=d, family=binomial))

summary(glm(y~width, data=d, family=binomial))

fit.col=glm(y~color, data=d, family=binomial)
anova(fit.null,fit.col,test="F")
anova(fit.null,fit.col,test="LRT")
summary(fit.col)

fit.spi=glm(y~spine, data=d, family=binomial)
anova(fit.null,fit.spi,test="F")
anova(fit.null,fit.spi,test="LRT")
summary(fit.spi)

fit.full=glm(y~weight+width+color+spine, data=d, family=binomial)

##stepwise BIC
res=step(fit.null, scope=list(lower=fit.null, upper=fit.full), direction="both", k=log(173))
##stepwise AIC
res=step(fit.null, scope=list(lower=fit.null, upper=fit.full), direction="both", k=2)

##Extract odds ratio
library(oddsratio)
or_glm(data=d, model=res, incr= list(width=1))

##goodness-of-fit test
library(LogisticDx)
g1=gof(res)
g1
unclass(g1)

library(ResourceSelection)
hoslem.test(res$y, fitted(res))

##Question 3
fit.null=glm(sat~1, data=d, family=poisson)

summary(glm(sat~weight, data=d, family=poisson))

summary(glm(sat~width, data=d, family=poisson))

fit.col=glm(sat~color, data=d, family=poisson)
anova(fit.null,fit.col,test="F")
anova(fit.null,fit.col,test="LRT")
anova(fit.null,fit.col,test="Chisq")
summary(fit.col)

fit.spi=glm(sat~spine, data=d, family=poisson)
anova(fit.null,fit.spi,test="F")
anova(fit.null,fit.spi,test="LRT")
summary(fit.spi)

fit.full=glm(sat~weight+width+color+spine, data=d, family=poisson)

##stepwise BIC
res=step(fit.null, scope=list(lower=fit.null, upper=fit.full), direction="both", k=log(173))
##stepwise AIC
res=step(fit.null, scope=list(lower=fit.null, upper=fit.full), direction="both", k=2)

with(res, cbind(res.deviance = deviance, df = df.residual,
  p = pchisq(deviance, df.residual, lower.tail=FALSE)))

library(MASS)
fit.nb=glm.nb(sat~weight+color, data=d)
fit.nb
summary(fit.nb)
library(pscl)
odTest(fit.nb)