#######################################################################
#  Revised nr.logit function that returns residuals, fitted values,
#  leverage and goodness-of-fit statistics
#######################################################################
nr.logit<-function(x,y,n,maxits=20,eps=1e-10,beta.start){
   if( missing(beta.start) ){
      emp.logit<-log((y+.5)/(n-y+.5))
      newbeta<-lsfit(x,emp.logit,intercept=F)$coef}
   else{newbeta<-beta.start}
   iter<-0
   converged<-F
   while( (!converged) & (iter<maxits) ){
      iter<-iter+1
      cat(iter); cat("...")
      beta<-newbeta
      tmp<-exp(x%*%beta)
      pi<-tmp/(1+tmp)
      mu<-n*pi
      w<-as.vector(n*pi*(1-pi))   # this is a vector, not a matrix
      xtwx<-t(w*x)%*%x
      xtwxinv<-solve(xtwx)
      newbeta<-beta + xtwxinv%*%t(x)%*%(y-mu)
      converged<-all(abs(newbeta-beta)<eps)}
   cat("\n")
   tmp<-exp(x%*%newbeta)
   pi<-as.vector(tmp/(1+tmp))
   loglik<-sum( y*log(pi) + (n-y)*log(1-pi) )
   # compute X2, Pearson residuals, influence
   r<-(y-n*pi)/sqrt(n*pi*(1-pi))
   X2<-sum(r^2)
   xstar<-diag(sqrt(pi*(1-pi)))%*%x
   hatstar<-xstar%*%solve(t(xstar)%*%xstar)%*%t(xstar)
   leverage<-diag(hatstar)
   # compute deviance contributions and G2
   deviance<-y*log(y/(n*pi)) + (n-y)*log((n-y)/(n-n*pi))
   if(any(y==0))
      deviance[y==0]<-n[y==0]*log(n[y==0]/(n[y==0]-n[y==0]*pi[y==0]))
   if(any(y==n))
      deviance[y==n]<-n[y==n]*log(n[y==n]/(n[y==n]*pi[y==n]))
   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=pi,
      expected.success=n*pi, expected.failure=n*(1-pi))
   result}

#######################################################################
#  Revised logit.print function that prints out goodness-of-fit 
#  information and allows a scale parameter
#######################################################################
logit.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")}
   tmp1<-sum(result$expected.success<5)
   tmp2<-sum(result$expected.failure<5)
   pct<-round(100*(tmp1+tmp2)/(2*length(y)),2)
   cat(paste(format(pct),
      "% of cells have expected counts below 5.0",sep=""))
   cat("\n")
   invisible()}

# Test the functions with the dose-response data
z<-c(-5,-4,-3,-2,-1)
x<-cbind(int=1,logdose=z)
y<-c(0,1,4,6,6)
n<-rep(6,5)
result<-nr.logit(x,y,n)
logit.print(result)