.packageName <- "glmmML"
glmmML <- function(formula,
                   family = binomial,
                   data,
                   cluster,
                   subset,
                   na.action,
                   offset,
                   start.coef = NULL,
                   start.sigma = NULL,
                   control = glm.control(epsilon = 1.e-8,
                       maxit = 100, trace = FALSE),
                   n.points = 16){

    method <- 1 ## Always vmmin! 1 if vmmin, 0 otherwise
    if (!method) stop("Use default method (the only available at present)")
    cl <- match.call()

    if (is.character(family)) 
        family <- get(family)
    if (is.function(family)) 
        family <- family()
    if (is.null(family$family)) {
        print(family)
        stop("`family' not recognized")
    }
    
    if (missing(data))
        data <- environment(formula)
    
    mf <- match.call(expand.dots = FALSE)
    ## get a copy of the call; result: a list.
    
    mf$family <- mf$start.coef <- mf$start.sigma <- NULL
    mf$control <- mf$maxit <- mf$boot <- NULL
    mf$n.points <- mf$method <- mf$start.coef <- mf$start.sigma <- NULL
    mf[[1]] <- as.name("model.frame") # turn into a call to model.frame
    mf <- eval(mf, environment(formula)) # run model.frame
    
    ## Pick out the parts.
    mt <-  attr(mf, "terms")
    
    
    xvars <- as.character(attr(mt, "variables"))[-1]
    if ((yvar <- attr(mt, "response")) > 0) 
        xvars <- xvars[-yvar]
    xlev <- if (length(xvars) > 0) {
        xlev <- lapply(mf[xvars], levels)
        xlev[!sapply(xlev, is.null)]
    }
    
    X <- if (!is.empty.model(mt)) 
        model.matrix(mt, mf, contrasts)
    
    p <- NCOL(X)
    
    Y <- model.response(mf, "numeric")
    offset <- model.offset(mf)
 
    cluster <- mf$"(cluster)"
    ##    return(clus)
    
    if (NCOL(Y) >  1) stop("Response must be univariate")
    
    if (!is.null(offset) && length(offset) != NROW(Y)) 
        stop(paste("Number of offsets is", length(offset), ", should equal", 
                   NROW(Y), "(number of observations)"))
    
    mixed <- ( !is.null(cluster) ) && ( n.points >= 2 )
    
    fit <- glmmML.fit(X, Y,
                      start.coef,
                      start.sigma,
                      mixed,
                      cluster,
                      offset,
                      family,
                      n.points,
                      control,
                      method,
                      intercept = ( attr(mt, "intercept") > 0)
                      )
    
    if (!fit$convergence) return(list(convergence = fit$convergence))
    bdim <- p + 1
    res <- list()

    res$boot <- FALSE
    res$convergence <- as.logical(fit$convergence)
    res$aic <- -2 * fit$loglik + 2 * (p + as.integer(mixed))
    res$variance <- fit$coef.variance
    if (mixed){
        res$sigma <- fit$sigma
        res$sigma.sd <- sqrt(fit$sigma.vari)
    }else{
        res$sigma = 0
        res$sigma.sd = 0
    }
    res$coefficients <- fit$beta
    res$deviance <- fit$deviance
    ##   options(show.error.messages = FALSE)
    ##  vari <- try(solve(-res$hessian))
    ##  if(is.numeric(vari)){
    ##    se <- sqrt(diag(vari))
    ##  }else{
    ##    se <- rep(NA, p + 1)
    ##  }
    res$df.residual <- fit$df.residual
    res$sd <- sqrt(diag(res$variance))
    names(res$sd) <- names(res$coefficients)
    res$mixed <- mixed
    if (mixed){
        res$frail <- fit$frail
    }
    res$call <- cl
    names(res$coefficients) <- c(colnames(X))
    class(res) <- "glmmML"
    res
}

glmmML.fit <- function (X, Y, 
                        start.coef = NULL, 
                        start.sigma = NULL,
                        mixed = FALSE,
                        cluster = NULL,                        
                        offset = rep(0, nobs),
                        family = binomial(),
                        n.points = 16,
                        control = glm.control(),
                        method,
                        intercept = TRUE,
                        boot = 0){
  
  X <- as.matrix(X)
  conv <- FALSE
  nobs <- NROW(Y)
  p <- NCOL(X)
  nvars <- p + as.integer(mixed)
  
  if (is.null(offset)) 
    offset <- rep(0, nobs)
  variance <- family$variance
  dev.resids <- family$dev.resids
  aic <- family$aic
  linkinv <- family$linkinv
  mu.eta <- family$mu.eta
  
  if (!is.function(variance) || !is.function(linkinv)) 
    stop("illegal `family' argument")

  if (is.null(start.coef)){
    start.coef <- numeric(p) # Start values equal to zero,
    if (family$family == "binomial"){
      start.coef[1] <- log(mean(Y) / (1 - mean(Y)))
    }else if (family$family == "poisson"){
      start.coef[1] <- log(mean(Y))
    }else{ ## this is a proviso!!
      start.coef[1] <- mean(Y)
    }
                           
  }else{                   
    if (length(start.coef) != p) stop("beta.start has wrong length")
  }
  
  if (mixed) {
    if (is.null(start.sigma)){
      start.sigma <- 0.5 ## More sofisticated choice is = ?
    }else{                  
      if (length(start.sigma) != 1) stop("sigma.start has wrong length")
    }
  }else{
    if (length(start.coef) != p) stop("beta.start has wrong length")
    n.points <- 1
    if (is.null(start.sigma)) start.sigma <- 0
  }
  
  ord <- order(cluster)
  Y <- Y[ord]
  X <- X[ord, ,drop = FALSE]

  ## Center the covariates so we avoid (some) numeric problems: 
  if (intercept){
    if (p >= 2){
      means <- numeric(p-1)
      for (i in 2:p){
        means[i-1] <- mean(X[, i])
        X[, i] <- X[, i] - means[i-1]
      }
    }
  }else{
    means <- numeric(p)
    for (i in 1:p){
      means[i] <- mean(X[, i])
      X[, i] <- X[, i] - means[i]
    }
  }
  
  cluster <- cluster[ord]
  fam.size <- as.vector(table(cluster))
  n.fam <- length(fam.size)
  

 if (family$family == "binomial"){
    if (family$link == "logit"){
      fam <- 0
    }else if (family$link == "cloglog"){
      fam <- 1
    }else{
      stop("Unknown link function; only 'logit' and 'cloglog' implemented")
    }
  }else if (family$family == "poisson"){
    fam <- 2
  }else{
    stop("Unknown family; only 'binomial' and 'poisson' implemented")
  }
              
  fit <- .C("glmm_ml",
            as.integer(fam),
            as.integer(method),
            as.integer(p),
            as.double(start.coef),
            as.double(start.sigma),
            as.double(t(X)),       ### Note CAREFULLY (03-01-09)!!!
            as.integer(Y),
            as.double(offset),
            as.integer(fam.size),
            as.integer(n.fam),
            as.integer(n.points),
            as.double(control$epsilon),
            as.integer(control$maxit),
            as.integer(control$trace),
            beta = double(p),  ## Return values from here.
            sigma = double(1),
            loglik = double(1),
            variance = double((p + 1) * (p + 1)),
            frail = double(n.fam),
            mu = double(nobs),
            convergence = integer(1),
            PACKAGE = "glmmML"
            )  

  vari <- matrix(fit$variance, ncol = (p + 1))
  ## Correct the estimate of the intercept for the centering:
  if (intercept){
    if (p >= 2){
      fit$beta[1] <- fit$beta[1] - sum(fit$beta[2:p] * means)
      aa <- numeric(p)
      aa[1] <- 1.0
      for (i in 2:p){ 
        aa[i] <- -means[i-1]
        ## "Restore" X (to what use?!):
        X[, i] <- X[, i] + means[i-1]
      }
      vari[1, 1] <- aa %*% vari[1:p, 1:p] %*% aa
    }
  }else{
    for (i in 1:p){ ## "Restore" mm (to what use?!):
      X[, i] <- X[, i] + means[i]
    }
  }
  
  if (mixed){
    sigma.vari <- vari[(p + 1), (p + 1)] * fit$sigma * fit$sigma
  }else{
    sigma.vari <- NULL
  }
  
  aic.model <- -2 * fit$loglik + 2 * nvars

  list(beta = fit$beta,
       sigma = fit$sigma,
       loglik = fit$loglik,
       coef.variance = vari[1:p, 1:p, drop = FALSE],
       sigma.variance = sigma.vari,
       frail = fit$frail,
       residuals = residuals,
       fitted.values = fit$mu, 
       family = family, 
       deviance = -2*fit$loglik,
       aic = aic.model, 
       #null.deviance = nulldev,
       df.residual = NROW(Y) - NCOL(X) - as.integer(mixed),
       df.null = NROW(Y) - as.integer(intercept),
       #y = y,
       convergence = fit$convergence)
}
glmmbootFit <- function (X, Y, 
                         start.coef = NULL, 
                         cluster = rep(1, length(Y)),                        
                         offset = rep(0, length(Y)),
                         family = binomial(),
                         control = glm.control(),
                         method,
                         boot, fortran = TRUE){

    X <- as.matrix(X)

    if (is.null(offset)) offset <- rep(0, length(Y))
    p <- ncol(X)
    if (is.null(start.coef)){
        start.coef <- numeric(p) # Start values equal to zero,
        if (FALSE){
        if (family$family == "binomial"){
            start.coef[1] <- log(mean(Y) / (1 - mean(Y)))
        }else if (family$family == "poisson"){
            start.coef[1] <- log(mean(Y))
        }else{ ## this is a proviso!!
            start.coef[1] <- mean(Y)
        }
    }
    }else{                   
        if (length(start.coef) != p) stop("beta.start has wrong length")
    }

    ord <- order(cluster)

    Y <- Y[ord]
    X <- X[ord, ,drop = FALSE]
    cluster <- cluster[ord]

    if (family$family == "binomial"){
        if (family$link == "logit"){
            fam <- 0
        }else if (family$link == "cloglog"){
            fam <- 1
        }else{
            stop("Unknown link function; only 'logit' and 'cloglog' implemented")
        }
    }else if (family$family == "poisson"){
        fam <- 2
    }else{
        stop("Unknown family; only 'binomial' and 'poisson' implemented")
    }
    
    famSize <- as.vector(table(cluster))
    nFam <- length(famSize)
  
    ## cat("nFam = ", nFam, "\n")

    if (fortran){
        if (NCOL(X)){
            means <- numeric(p)
            means <- colMeans(X)
            X <- scale(X, center = TRUE, scale = FALSE)
    
            fit <- .C("glmm_boot",
                      as.integer(fam),
                      as.integer(method),
                      as.integer(p),
                      as.double(start.coef),
                      as.integer(cluster),
                      as.double(t(X)),       ## Note! ##
                      as.integer(Y),
                      as.double(offset),
                      as.integer(famSize),
                      as.integer(nFam),
                      as.double(control$epsilon),
                      as.integer(control$maxit),
                      as.integer(control$trace),
                      as.integer(boot),
                      beta = double(p),
                      loglik = double(1),
                      hessian = double(p * p),
                      frail = double(nFam),
                      bootP = double(1),
                      bootLog = double(boot),
                      convergence = integer(1)
                      )
            res <- list(coefficients = fit$beta,
                        logLik = fit$loglik,
                        frail = fit$frail,
                        bootLog = fit$bootLog,
                        bootP = fit$bootP)
            res$variance <- solve(-matrix(fit$hessian, nrow = p, ncol = p))
            res$sd <- sqrt(diag(res$variance))
            res$boot_rep <- boot

            return(res)
        }else{ # A null model:
            ## Center the covariates so we avoid (some) numeric problems: 
            fit <- .C("glmm_boot0",
                      as.integer(fam),
                      as.integer(method),
                      ##as.integer(p),
                      ##as.double(start.coef),
                      as.integer(cluster),
                      ##as.double(t(X)),       ## Note! ##
                      as.integer(Y),
                      as.double(offset),
                      as.integer(famSize),
                      as.integer(nFam),
                      ##as.double(control$epsilon),
                      ##as.integer(control$maxit),
                      as.integer(control$trace),
                      as.integer(boot),
                      ##beta = double(p),
                      loglik = double(1),
                      ##hessian = double(p * p),
                      frail = double(nFam),
                      bootP = double(1),
                      bootLog = double(boot),
                      convergence = integer(1)
                      )
            res <- list(coefficients = NULL,
                        logLik = fit$loglik,
                        frail = fit$frail,
                        bootLog = fit$bootLog,
                        bootP = fit$bootP)
            res$variance <- NULL
            res$sd <- NULL
            res$boot_rep <- boot

            return(res)
        }

    }else{
        
        fit <- snut(Y, X, cluster, offset)
        
        if (boot < 1) error("boot must be at least 1")
        logl <- numeric(boot + 1)
        logl[1] <- -fit$value
        
        for (i in 2:(boot + 1)){
            cat("****************** Replicate ", i-1, "\n")
            cluster <- sample(cluster)
            ord <- order(cluster)
            cluster <- cluster[ord]
            Y <- Y[ord]
            X <- X[ord, ,drop = FALSE]
            
            logl[i] <- -snut(Y, X, cluster, offset)$value
        }
        
        ##cat("logl = ", logl, "\n")
        p.value <- 1 - rank(logl, ties.method = "first")[1] / (boot + 1)  
        
        nvars <- length(unique(cluster))
        nvar <- length(fit$par)
        aic.model <- 2 * fit$value + 2 * nvars
        ## Note difference between nvar & nvars!
        n.fam <- nvar - p
        list(beta = fit$par[(n.fam + 1):nvar],
             loglik = -fit$value,
             ##coef.variance = vari[1:p, 1:p, drop = FALSE],
             frail = fit$par[1:n.fam],
             ##residuals = residuals,
             ##fitted.values = fit$mu, 
             ##family = family, 
             deviance = 2*fit$value,
             aic = aic.model, 
             ##null.deviance = nulldev,
             ##df.residual = NROW(Y) - NCOL(X) - as.integer(mixed),
             ##df.null = NROW(Y) - as.integer(intercept),
             ##y = y,
             ##convergence = fit$convergence
             frail.p = p.value
             )
    }
}
glmmboot <- function(formula,
                   family = binomial,
                   data,
                   cluster,
                   subset,
                   na.action,
                   offset,
                   start.coef = NULL,
                   control = glm.control(epsilon = 1.e-8,
                       maxit = 100, trace = FALSE),
                   boot = 0){

    method <- 1 ## Always vmmin! 1 if vmmin, 0 otherwise
    if (!method) stop("Use default method (the only available at present)")
    cl <- match.call()

    if (is.character(family)) 
        family <- get(family)
    if (is.function(family)) 
        family <- family()
    if (is.null(family$family)) {
        print(family)
        stop("`family' not recognized")
    }
    
    if (missing(data))
        data <- environment(formula)
    
    mf <- match.call(expand.dots = FALSE)
    ## get a copy of the call; result: a list.
    
    mf$family <- mf$start.coef <- mf$start.sigma <- NULL
    mf$control <- mf$maxit <- mf$boot <- NULL
    mf$n.points <- mf$method <- mf$start.coef <- NULL
    mf[[1]] <- as.name("model.frame") # turn into a call to model.frame
    mf <- eval(mf, environment(formula)) # run model.frame
    
    ## Pick out the parts.
    mt <-  attr(mf, "terms")
    
    
    xvars <- as.character(attr(mt, "variables"))[-1]
    if ((yvar <- attr(mt, "response")) > 0) 
        xvars <- xvars[-yvar]
    xlev <- if (length(xvars) > 0) {
        xlev <- lapply(mf[xvars], levels)
        xlev[!sapply(xlev, is.null)]
    }
    
    X <- if (!is.empty.model(mt)) 
        model.matrix(mt, mf, contrasts)
    
    p <- NCOL(X)
    
    Y <- model.response(mf, "numeric")
    offset <- model.offset(mf)
 
    cluster <- mf$"(cluster)"
    
    if (NCOL(Y) >  1) stop("Response must be univariate")
    
    if (!is.null(offset) && length(offset) != NROW(Y)) 
        stop(paste("Number of offsets is", length(offset), ", should equal", 
                   NROW(Y), "(number of observations)"))
    
    ## Remove eventual intercept from X.
    ## Taken care of thru separate intercepts for each 'cluster'.
    if (!is.na(coli <- match("(Intercept)", colnames(X))))
        X <- X[, -coli, drop = FALSE]

    fortran <- TRUE
    res <- glmmbootFit(X, Y,
                       start.coef,
                       cluster,
                       offset,
                       family,
                       control,
                       method,
                       boot,
                       fortran)
    
    res$mixed <- FALSE
    res$deviance <- -2 * res$logLik
    nvars <- NCOL(X) + length(unique(cluster))
    res$df.residual <- length(Y) - nvars
    res$aic <- res$deviance + 2 * nvars
    res$boot <- TRUE
    res$call <- cl
    names(res$coefficients) <- c(colnames(X))
    class(res) <- "glmmboot"
    res
}

print.glmmML <- function(x,
                         digits = max(3, getOption("digits") - 3),
                         na.print = "",
                         ...){ 
    
    cat("\nCall: ", deparse(x$call), "\n\n")
    savedig <- options(digits = digits)
    on.exit(options(savedig))
    coef <- x$coefficients
    se <- x$sd
    tmp <- cbind(coef,
                 se,
                 coef/se,
                 signif(1 - pchisq((coef/se)^2, 1), digits - 1)
                 )
    dimnames(tmp) <- list(names(coef),
                          c("coef", "se(coef)", "z", "Pr(>|z|)")
                          )
    cat("\n")
    prmatrix(tmp)
    
    if(x$mixed){
        cat("\nStandard deviation in mixing distribution: ", x$sigma,  "\n")
        cat("Std. Error:                                ", x$sigma.sd, "\n")
    }
    if(x$boot){
        cat("\n Bootstrap p-value for fixed mixing: ",
        x$bootP, "(", x$boot_rep, ")\n")
    }
    cat("\nResidual deviance:",
        format(signif(x$deviance, digits)), "on",
        x$df.residual, "degrees of freedom", 
        "\tAIC:",
        format(signif(x$aic, digits)), "\n")
}
print.glmmboot <- function(x,
                         digits = max(3, getOption("digits") - 3),
                         na.print = "",
                         ...){ 
    
    cat("\nCall: ", deparse(x$call), "\n\n")
    savedig <- options(digits = digits)
    on.exit(options(savedig))

    if (length(x$coefficients)){
        coef <- x$coefficients
        se <- x$sd
        tmp <- cbind(coef,
                     se,
                     coef/se,
                     signif(1 - pchisq((coef/se)^2, 1), digits - 1)
                     )
        dimnames(tmp) <- list(names(coef),
                              c("coef", "se(coef)", "z", "Pr(>|z|)")
                              )
        cat("\n")
        prmatrix(tmp)
    }
    
    if (x$boot_rep){
        cat("\n Bootstrap p-value for fixed mixing: ",
            x$bootP, "(", x$boot_rep, ")\n")
    }
    
    cat("\nResidual deviance:",
        format(signif(x$deviance, digits)), "on",
        x$df.residual, "degrees of freedom", 
        "\tAIC:",
        format(signif(x$aic, digits)), "\n")
}
snut <- function(Y, X,
                 cluster = rep(1, length(Y)),
                 offset = rep(0, length(Y))
                 ){
    ## X may not have a constant column!
    n <- length(Y)
    if (is.matrix(X)){
        if (NROW(X) != n) stop("Wrong dimension of X")
        q <- NCOL(X)
    }else{
        if (length(X) != n) stop("Wrong length of X")
        q <- 1
        X <- matrix(X, ncol = 1)
    }

    cluster <- as.vector(unclass(factor(cluster)))

    in.here <- function(y) (sum(y) > 0.5) && (sum(y) < (length(y) - 0.5))
    
    here <- cluster %in% which(tapply(Y, cluster, in.here))
    Y <- Y[here]
    X <- X[here, , drop = FALSE]
    cluster <- cluster[here]
    offset <- offset[here]
    
    p <- length(unique(cluster))
    pq <- p + q
    beta <- numeric(pq)

    cluster <- as.vector(unclass(factor(cluster)))
    
    fun <- function(beta){
        ## Minus The log likelihood:
        lin <- offset + beta[cluster] + X %*% beta[(p+1):pq]
        -sum(Y * lin - log(1 + exp(lin)))
    }

    gra <- function(beta){
        # Minus the gradient
        ret <- numeric(length(beta))
        lin <- offset + beta[cluster] + X %*% beta[(p+1):pq]
        P <- exp(lin)
        P <- P / (1 + P)
        ymP <- Y - P # is an nx1 matrix
        ret[1:p] <- tapply(ymP, cluster, sum)
        ret[(p+1):pq] <- t(ymP) %*% X
        ##for (s in 1:q) ret[s + p] <- sum(X[, s] * ymP)
        -ret
    }

    res <- optim(beta, fun, gra, method = "BFGS", control = list(trace = TRUE))
    res
}
summary.glmmML <- function(object, ...){
    print.glmmML(object, ...)
}
summary.glmmboot <- function(object, ...){
    print.glmmboot(object, ...)
}
