######################################################################
#  Arguments
#           x = vector or matrix of covariates
#           y = vector of responses ("number of successes")
#           n = vector of binomial indices ("number of trials")
#      maxits = maximum number of iterations (default is 20)
#         eps = convergence criterion (default is 1E-10)
#  beta.start = optional starting value for beta
#
#  Value returned
#     A list with the following elements:
#        beta = final estimate
#    cov.beta = estimated covariance matrix
#        iter = number of iterations
#   converged = T if converged, F otherwise
#      loglik = final value of loglikelihood
#
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<-tmp/(1+tmp)
   loglik<-sum( y*log(pi) + (n-y)*log(1-pi) )
   result<-list(beta=newbeta, cov.beta=xtwxinv, iter=iter, 
      converged=converged, loglik=loglik)
   result}

#######################################################################
#  Function for nicely formatted output
#######################################################################
logit.print<-function(result){
   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(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")
   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)