ML-NR

# author andrew.marete@canada.ca
# Generalized linear model: Maximum likelihood with Newton-Raphson

ML_NR <- function(y=NULL, x=NULL, th=1e-6, start_theta=c(0,0), maxIt=100){
  # y: vector of  observations
  # x: vector of covariates
  # th: convergence threshold, i.e., iterative process stops
  #     when delta theta < th
  # start_theta: starting parameter
  # maxIt : The maximum number of iterations

  theta_t <- start_theta
  it <- d.mu <- d.b <- 1
  theta_out <- NULL
  while(( abs(d.b) > th | abs(d.mu) > th ) & it < maxIt) {
    b <- theta_t[1]
    mu <- theta_t[2]

    S <- matrix(ncol=1, nrow=2)
    H <- matrix(ncol=2, nrow=2)
    rownames(S) <- rownames(H) <- colnames(H) <- c("b", "mu")

    Logit <- exp(b*x+mu)/(1+exp(b*x+mu))
    Logit2 <- exp(b*x+mu)/(1+exp(b*x+mu))^2

    S["mu",] <- sum(y - Logit)
    S["b",] <- sum(x*y - Logit*x)

    H["mu","mu"] <- -sum(Logit2)
    H["mu", "b"] <- -sum(Logit2*x)
    H["b","b"] <- -sum(Logit2*x**2)
    H["b", "mu"] <- -sum(Logit2*x)

    iH <- solve(H)
    theta_u <- theta_t - iH%*%S
    d.b <- theta_u[1] - theta_t[1]
    d.mu <- theta_u[2] - theta_t[2]

    theta_t <- theta_u

    theta_out <- rbind(theta_out, c(t(theta_u),d.b, d.mu))
    colnames(theta_out) <- c("beta", "mu", "delta_beta", "delta_mu")
    it <- nrow(theta_out)

  }
  theta <- theta_out
  vc_theta <- -iH

  # - - - #
  layout(matrix(c(1,2), ncol=1, byrow=T), heights=c(5,.4))
  par(mar=c(5,4,4,3))
  plot(x=1:nrow(theta), y=theta[,3], bty="n", las=1, type="l", lwd=2, lty=5, ylab=expression(paste(Delta, "(", beta,")")), xlab="Iterations",main = "Convergence Criteria - Newton Raphson", cex.axis=.8)
  par(new=TRUE)
  plot(x=1:nrow(theta), y=theta[,4], bty="n", las=1, type="l", lwd=2, lty=3, yaxt="n", xlab="", ylab="", xaxt="n")
  axis(side=4, cex.axis=.8, las=1)
  mtext(side=4, line=2, expression(paste(Delta, "(", mu,")")))

  par(mai=c(0,0,0,0))
  plot.new()

  legend("center", legend=c(expression(paste(Delta, "(", mu,")")),
                            expression(paste(mu, "(", beta,")"))), lty=c(5,3), cex=.8, bty="n", horiz=TRUE)
  # - - - #

  res <- list(theta=theta, "Assy variance-covariance Matrix"=vc_theta)
  return(res)
}

file = "https://www.dropbox.com/s/f6pt59jtnrfzufi/dataAlbert2Col.txt?dl=1"
dat <- read.delim(file, sep=" ")
fit <- ML_NR(y = dat[,1], x = dat[,2])
fit