# Edps 589
# Fall 2018
# C.J.Andersion
#
#  Answer for Homework 5, problem 2 (3.13 & 3.14 from Agresti)
#
library(MASS)

setwd("C:\\Users\\cja\\Dropbox\\edps 589\\homework")

crabs<- read.table("crab_data.txt",header=TRUE)

##############################################
# Mean and variance of number of satellites  #
############################################## 
mean(crabs$satell)
var(crabs$satell)
length(crabs$satell)
summary(crabs$satell)

##################################################
# 3.13 from Agresti                              #
##################################################

##############################################
# Histogram of number of satellites          #
##############################################

hist(crabs$satell,
breaks=16,
xlab="Number of Satellites",
main="Distribution of Number of Satellites",col="skyblue"
)
text(12,80,"mean = 2.92")
text(12,77,"  var = 0.01")
text(12,74,"   N = 173")
text(12,71,"  min = 0")
text(12,68,"median = 2")
text(12,65,"  max =  14")


plot(crabs$weight,jitter(crabs$satell),
     type="p",
	 pch=20,
	 main="Number of Satellites by Weight...a little jittering",
	 xlab="Weight of Crab",
	 ylab="Number of Satellites")
smooth <- loess(satell ~ weight,data=crabs)
j <- order(crabs$weight)
lines(crabs$weight[j],smooth$fitted[j],col="blue")	 


crabs$logweight <- log(crabs$weight)
plot(crabs$logweight,jitter(crabs$satell),
     type="p",
	 pch=20,
	 main="Number of Satellites by log(Weight)...a little jittering",
	 xlab="log(Weight of Crab)",
	 ylab="Number of Satellites")
smooth <- loess(satell ~ logweight,data=crabs)
j <- order(crabs$weight)
lines(crabs$logweight[j],smooth$fitted[j],col="blue")	


# First rescale to kg
crabs$weight <- crabs$weight/1000

# 1. Fit poisson log-linear model
summary( mod1 <- glm(satell ~ weight,data=crabs,family=poisson) )

# 2. Estimated mean if weight=2.44
(muhat_2.44<- exp(mod1$coefficients[1] + mod1$coefficients[2]*2.44))

# 3. For 1 unit change
exp(mod1$coefficients[2])

#  CI for beta using estimates
(lowerB <-  mod1$coefficients[2] - 1.96*0.06502)
(upperB <-  mod1$coefficients[2] + 1.96*0.06502)

#   CI for mean effect
exp(lowerB)
exp(upperB)

# 4. Wald
# either report z or 
9.064**2 = 82.16

1 - pchisq(82.16,1)

# 5. Likelihood ratio test

anova(mod1)

1-pchisq(71.925,1)


##################################################
# 3.14 from Agresti                              #
##################################################

summary( mod.nb <- glm.nb(satell ~ weight, data=crabs))

#sas reports 1/phi, which for these data is
1/0.9310

#  CI for beta using estimates
(lowerB <-  mod1$coefficients[2] - 1.96*0.06502)
(upperB <-  mod1$coefficients[2] + 1.96*0.06502)

#   CI for mean effect
exp(lowerB)
exp(upperB)

# 4. Wald
# either report z or 
4.817**2 = 23.20

1 - pchisq(23.20,1)

# 5. Likelihood ratio test

anova(mod.nb)

1-pchisq(20.27,1)



##################################################
# Zero inflated model                            #
##################################################

library(ggplot2)
library(pscl)

summary(zip <- zeroinfl(satell ~ weight | width , data = crabs))