logitreg <- function(y, x, start, tol = 10^-8, Nmax = 20) {

# *************************************************************************
#
# FILE: logitreg.r
#
# Date Created: Feb 24, 2004
# Author: Mark E. Irwin
#
# Contents:  Performs logistic regression by Fisher scoring.  Assumes that
#            the response variable is binary (0/1).  Automatically adds
#            a vector of ones to deal with the intercept term.  Instead of
#            using sums as described by Lange, matrix forms are used instead.
#            The error tolerance is for relative error
#
# Revision History
# Date          Name            Changes/Reasons
#
# *************************************************************************
#
# Input:
#
#   y: response vector
#   x: covariates
#   start: starting estimate of theta
#   tol: convergence tolerance (maximum relative error)
#   Nmax: maximum number of iterations
#
# Output:
#
#   List containing:
#
#     theta: estimate of regression parameters
#     info:  expected information matrix for parameter estimates
#     niter: number of Fisher scoring steps
#
# *************************************************************************

# Setup variables for scoring algorithm

  if(is.matrix(x))
    obs <- dim(x)[1]
  else
    obs <- length(x)

  ntheta <- length(start)

  Z <- cbind(rep(1,obs),x)
  thetaold <-matrix(start,ncol=1)
  thetanew <- thetaold

  n <- 0
  error <- 2 * tol

# Run loop

  while (error > tol && n < Nmax) {

    n <- n + 1
    thetaold <- thetanew
    
    logit <- Z %*% thetaold
    mu <- plogis(logit)
    qprime <- matrix(rep(dlogis(logit), ntheta), ncol=ntheta)
    dmu <- qprime * Z
    sigma2 <- mu * (1 - mu)
    sc <- matrix(rep((y - mu) / sigma2, ntheta), ncol=ntheta) * dmu
    score <- matrix(apply(sc, 2, sum), ncol = 1)

    vr <- dmu / matrix(rep(sqrt(sigma2), ntheta), ncol=ntheta)
    J <- t(vr) %*% vr

    thetanew <- thetaold + solve(J, score)
    err <- abs(thetanew - thetaold) / (abs(thetaold) + 0.1)
    error <- max(err)

  }

if (n > Nmax && error > tol)
  warning('Algorithm did not converge in desired number of steps.');

list(theta=thetanew, info=J, niter=n)
}