# Poisson regression example using the Newton-Raphson method
# Data Table 4.5 Number of cases of AIDS in Australia for
# successive quarters 1984-1988
y=c(1,6,16,23,27,39,31,30,43,51,63,70,88,97,91,104,110,113,149,159)
# For larger data sets use read.table function!
ind=c(1:20)
plot(ind,y,xlab="time",ylab="counts")
lines(ind,ind,col="green")
lines(ind,ind^1.5, col="blue")
lines(ind,ind^1.7,col="red")
lines(ind,ind^2,col="magenta")
plot(log(ind),log(y),xlab="log(time)",ylab="log(counts)")
abline(a=0,b=1.7,col="blue")
# Define elements for N-R recursions
#Model without intercept
n=length(y)
X=cbind(log(c(1:20)))
theta=3
iter=10
for(i in 1:iter)
{
w=exp(theta*log(c(1:20)) )
z=theta*log(c(1:20)) + (y - w )/w
W=diag(w)
tn=solve(t(X)%*%W%*%X)%*%t(X)%*%W%*%z
theta=tn
}

# Model with intercept
n=length(y)
X=cbind(rep(1,n), log(c(1:20)) )
b=c(1,1)
iter=10
for(i in 1:iter)
{
w=exp(b[1]+b[2]*log(c(1:20)) )
z=b[1]+ b[2]*log(c(1:20)) + (y - w )/w 
W=diag(w)
bn=solve(t(X)%*%W%*%X)%*%t(X)%*%W%*%z
b=bn
}
# Fit models directly with glm function
fit1=glm(y~log(ind),family=poisson(link="log"))
fit2=glm(y~log(ind)-1,family=poisson(link="log"))
# Try fit=glm(y~log(ind),family=poisson(link="identity"))