#########################
#
# File: 20111101.R
#       Generalized Linear Models (Poisson)
#       Count Dependent Variables
#
#########################


# Preamble

library(RFS)



# Let's read in some pseudo data

prv <- read.csv("http://courses.kvasaheim.com/pols6123/data/fakepoisson.csv")
names(prv)
attach(prv)
plot(x,y)


m1 <- glm(y~x, family=gaussian(link=log))


m1 <- glm(y~x, family=gaussian(link=log), start=c(0,0.5) )
summary(m1)


newX <- seq(-10, 40, length=100)
pr   <- exp( predict(m1, newdata=data.frame(x=newX) )  )
lines(newX,pr, col="tomato",lwd=3)

res  <- residuals(m1)
plot(x,res)
plot(y,res)

mean(res)



m2 <- glm(y~x, family=poisson(link=log))
summary(m2)

m3 <- glm(y~x, family=quasipoisson(link=log))
summary(m3)

newX <- seq(-10, 40, length=100)
pr   <- exp( predict(m3, newdata=data.frame(x=newX) )  )
lines(newX,pr, col="darkgreen",lwd=3)


##################################################


detach(prv)
rm(list=ls())

ter <- read.csv("http://courses.kvasaheim.com/pols6123/data/terrorism.csv")
names(ter)
summary(ter)
attach(ter)

# Model PIRA killings

m1p <- glm(pira ~ labour*riteleft, family=poisson(link=log))
summary(m1p)

m1q <- glm(pira ~ labour*riteleft, family=quasipoisson(link=log))
summary(m1q)

library(MASS)
m1n <- glm.nb(pira ~ labour*riteleft, link=log)
summary(m1n)

plot(riteleft,pira)

newX <- seq(-50,50,length=1000)
pr1 <- exp( predict(m1q, newdata=data.frame(riteleft=newX, labour=1) ) )
pr0 <- exp( predict(m1q, newdata=data.frame(riteleft=newX, labour=0) ) )

lines(newX,pr1, col=2)
lines(newX,pr0, col=4)




# Negative binomial model predictions

plot(riteleft,pira)

newX <- seq(-50,50,length=1000)
pr1 <- exp( predict(m1n, newdata=data.frame(riteleft=newX, labour=1) ) )
pr0 <- exp( predict(m1n, newdata=data.frame(riteleft=newX, labour=0) ) )

lines(newX,pr1, col=2)
lines(newX,pr0, col=4)



# New set of models

m2q <- glm(pira ~ labour*(riteleft+I(riteleft^2)), family=quasipoisson(link=log))
m2q <- glm(pira ~ labour*riteleft + I(riteleft^2), family=quasipoisson(link=log))
m2q <- glm(pira ~ labour + riteleft + I(riteleft^2), family=quasipoisson(link=log))
summary(m2q)

m2n <- glm.nb(pira ~ labour*(riteleft+I(riteleft^2)), link=log)
m2n <- glm.nb(pira ~ labour*riteleft + I(riteleft^2), link=log)
summary(m2n)


# Compare predictions

plot(riteleft,pira)
newX <- seq(-50,50,length=1000)

pr1 <- exp( predict(m2q, newdata=data.frame(riteleft=newX, labour=1) ) )
pr0 <- exp( predict(m2q, newdata=data.frame(riteleft=newX, labour=0) ) )
lines(newX,pr1, col="darkblue")
lines(newX,pr0, col="lightblue")

pr1 <- exp( predict(m2n, newdata=data.frame(riteleft=newX, labour=1) ) )
pr0 <- exp( predict(m2n, newdata=data.frame(riteleft=newX, labour=0) ) )
lines(newX,pr1, col="darkred")
lines(newX,pr0, col="pink")






m3q <- glm(pira ~ labour + riteleft+I(riteleft^2)+I(riteleft^3)+I(riteleft^4) , family=quasipoisson(link=log))
summary(m3q)

plot(riteleft,pira, xlim=c(-100,100))
newX <- seq(-50,50,length=1000)
pr1 <- exp( predict(m3q, newdata=data.frame(riteleft=newX, labour=1) ) )
pr0 <- exp( predict(m3q, newdata=data.frame(riteleft=newX, labour=0) ) )
lines(newX,pr1, col="darkblue",  lwd=2)
lines(newX,pr0, col="darkgreen", lwd=2)




############################################################


      Ten-minute break


############################################################

rm(list=ls())


set.seed(557)
time <- numeric()
dose <- round(runif(100)*9)+1
for(i in 1:100) time[i] <- round(rexp(1,1/dose[i]))+1

plot(dose,time, las=1, xlim=c(0,13), ylim=c(0,30))



m1 <- glm(time ~ dose, family=Gamma(link=inverse))
summary(m1)

newX <- seq(1,10, length=100)
pr <- 1/predict(m1, newdata=data.frame(dose=newX) )
lines(newX,pr , col=2, lwd=4)