#######################################################################
poisson.regression<-function(x,y,offset=0,maxits=20,eps=1e-10,beta.start){
   # Fit a Poisson model with a log link
   if( missing(beta.start) ){
      ystar<-log(y+.5)-offset
      newbeta<-lsfit(x,ystar,intercept=F)$coef}
   else{newbeta<-beta.start}
   if(is.vector(x)) x<-matrix(x,ncol=1)
   iter<-0
   converged<-F
   while( (!converged) & (iter<maxits) ){
      iter<-iter+1
      cat(iter); cat("...")
      beta<-newbeta
      eta<-as.vector(x%*%beta)+offset
      mu<-exp(eta)
      z<-eta - offset + (y - mu)/mu
      weights<-mu
      tmp<-lsfit(x,z,wt=weights,intercept=F)
      newbeta<-tmp$coef
      converged<-all(abs(newbeta-beta)<eps)}
   cat("\n")
   eta<-as.vector(x%*%beta)
   mu<-exp(eta)
   loglik<-sum( y*log(mu) - mu)
   xtwxinv<-ls.diag(tmp)$cov.unscaled
   # compute X2, Pearson residuals, influence
   r<-(y-mu)/sqrt(mu)
   X2<-sum(r^2)
   xstar<-diag(sqrt(mu))%*%x
   hatstar<-xstar%*%solve(t(xstar)%*%xstar)%*%t(xstar)
   leverage<-diag(hatstar)
   # compute deviance contributions and G2
   deviance<-y*log(y/mu)-(y-mu)
   if(any(y==0))
      deviance[y==0]<-mu[y==0]
   G2<-2*sum(deviance)
   result<-list(beta=newbeta, cov.beta=xtwxinv, iter=iter, 
      converged=converged, loglik=loglik, residuals=r,
      X2=X2, G2=G2, df=nrow(x)-ncol(x), leverage=leverage, fitted=mu,
      eta=eta, z=z)
   result}
#######################################################################
poisson.print<-function(result, scale=1.0){
   if(result$converged){
      cat(paste("The Newton-Raphson algorithm converged in",
         format(result$iter),"iterations."))}
   else{
      cat(paste("Algorithm failed to converge by",
         format(result$iter),"iterations"))}
   cat("\n\n")
   coef<-result$beta
   SE<-sqrt(scale*diag(result$cov.beta))
   tstat<-coef/SE
   pval<-2*pnorm(-abs(tstat))
   coef.table<-cbind(coef,SE,round(tstat,2),round(pval,3))
   dimnames(coef.table)<-list(names(result$beta),
       c("coef","SE","coef/SE","pval"))
   print(coef.table)
   cat("\n")
   cat(paste("Loglikelihood =",result$loglik))
   cat("\n")
   cat(paste("Pearson's X^2 =",result$X2))
   cat("\n")
   cat(paste(" Deviance G^2 =",result$G2))
   cat("\n")
   cat(paste("           df =",format(result$df)))
   cat("\n")
   if(scale!=1.0){
      cat(paste("   Scaled X^2 =",result$X2/scale))
      cat("\n")
      cat(paste("   Scaled G^2 =",result$G2/scale))
      cat("\n")}
   tmp<-sum(result$fitted<5)
   pct<-round(100*tmp/length(y),2)
   cat(paste(format(pct),
      "% of observations have expected counts below 5.0",sep=""))
   cat("\n")
   tmp<-min(result$fitted)
   cat(paste("The minimum expected cell count is ",
      format(tmp),sep=""))
   cat("\n")
   invisible()}
#######################################################################