R/ml_probit.R

In ui: Uncertainty Intervals and Sensitivity Analysis for Missing Data

#' Fit maximum likelihood for fixed values of rho
#'
#'This is a support function for \code{\link{ui_probit}}
#' @param out_formula Formula for outcome regression.
#' @param mis_formula Formula for regression model for the missingness mechanism.
#' @param data Data frame containing the variables in the formulas
#' @param rho Vector containing the values of rho for which we want to fit the likelihood.
#' @param progress If TRUE prints out process time for each maximazation of the likelihood.
#' @param method Maximazation method to be passed through \code{maxLik}
#' @importFrom maxLik maxLik
#' @importFrom stats get_all_vars glm binomial update.formula vcov setNames as.formula
#' @export
ml_probit<-function(out_formula,mis_formula,data,rho,progress=TRUE,method='NR'){

  ###Warnings
if(!inherits(data,"data.frame")){
	stop('Data must be a data frame')
    }

  y<-get_all_vars(out_formula,data=data)[,1]
  z <- as.numeric(!is.na(y))

  data$z<-z
  out_model<-glm(out_formula,family=binomial(link="probit"),data=data)

  m_f<-update.formula(out_formula, as.formula(paste('z', "~ .")))
	#Need to keep all NA in Xy, if different covariates in mis_model and out_model 	make new model without missing in y
  if(is.null(mis_formula)){
  	mis_formula<-m_f
	mis_model<-glm(mis_formula,family=binomial(link="probit"),data=data)
	t_m<-mis_model
    }else{
	mis_model<-glm(mis_formula,family=binomial(link="probit"),data=data)
	t_m<-glm(m_f,family=binomial(link="probit"),data=data)
}

  Xy<-model.matrix(t_m)
  Xz<-model.matrix(mis_model)

  nrho <- length(rho)
  p1 <- length(coef(out_model))
  p2 <- length(coef(mis_model))
  d<- p1+p2

  coef <- matrix(nrow=(p1 + p2), ncol = nrho)
  rownames(coef) <- c(names(out_model$coef), names(mis_model$coef))
  colnames(coef)<-paste(rho)

  value <- vector(length = nrho)
  names(value) <- paste(rho)
  max_type <- value

  i0 <- which(rho == 0)
  coef[, i0] <- c(coef(out_model), coef(mis_model))

  max.info <- list(converged = array(dim=nrho), message = array(dim=nrho), method = array(dim=nrho),
                   iterations = array(dim=nrho))
  max.info$converged[i0] <- mis_model$converged == TRUE & out_model$converged == TRUE
  max.info$message[i0] <- " "
  max.info$method[i0] <- "glm"
  max.info$iterations[i0] <- NA


  glmVcov <- vcov(out_model)
  mis_glmVcov <- vcov(mis_model)

  value[i0] <- logl_probit(coef[, i0], rho = 0, Xz = Xz, Xy = Xy, y = y, z = z)

  sigma <- list()
  sigma[[i0]] <- matrix(0, nrow = d, ncol = d)
  sigma[[i0]][1:p1, 1:p1] <- glmVcov
  sigma[[i0]][(p1 + 1):d, (p1 + 1):d] <- mis_glmVcov
  sigma[[i0]][(p1 + 1), (p1 + 1)] <- summary(out_model)$dispersion^2/(2*out_model$df.residual)
  dimnames_sigma <- c(dimnames(glmVcov)[[1]], dimnames(mis_glmVcov)[[1]])
  dimnames(sigma[[i0]]) <- list(dimnames_sigma, dimnames_sigma)

  if(sum(rho < 0) > 0){
    for(i in 1:(i0 - 1)){

      if(progress == TRUE){
        cat("Optimization for rho =", rho[i0 - i], sep=" ", "\n")
        ptm <- proc.time()
      }

      f <- function(par)
        logl_probit(par, rho = rho[i0 - i], Xz = Xz, Xy = Xy, y = y, z=z)

      g <- function(par)
        grr(par, rho = rho[i0 - i], Xz = Xz, Xy = Xy, y = y, z=z)

      h <- function(par)
        hess(par, rho = rho[i0 - i], Xz = Xz, Xy = Xy, y = y, z=z)

      ML <- maxLik(logLik = f, grad = g, hess = h, start = coef[, i0 - i + 1], method = method)

      max.info$converged[i0 - i] <- ML$code <= 2
      max.info$message[i0 - i] <- ML$message
      max.info$method[i0 - i] <- ML$type
      max.info$iterations[i0 - i] <- ML$iterations
      coef[,i0 - i] <- ML$estimate
      value[i0 - i] <- ML$maximum

      sigma[[i0 - i]] <- -solve(ML$hessian)


      if(progress == TRUE)
        cat("   Time elapsed:", (proc.time()-ptm)[3], "s", "\n")
    }
  }

  if(sum(rho > 0) > 0){
    for(i in (i0 + 1):nrho){

      if(progress == TRUE){
        cat("Optimization for rho =", rho[i], sep=" ", "\n")
        ptm <- proc.time()
      }

      f <- function(par)
        logl_probit(par, rho = rho[i], Xz = Xz, Xy = Xy, y = y, z = z)

      g <- function(par)
        grr(par, rho = rho[i], Xz = Xz, Xy = Xy, y = y , z = z)

      h <- function(par)
        hess(par, rho = rho[i], Xz = Xz, Xy = Xy, y = y, z = z)

      ML <- maxLik(logLik = f, grad = g, hess = h, start = coef[, i - 1], method = method)

      max.info$converged[i] <- ML$code <= 2
      max.info$message[i] <- ML$message
      max.info$method[i] <- ML$type
      max.info$iterations[i] <- ML$iterations
      coef[,i] <- ML$estimate
      value[i]  <-  ML$maximum

      sigma[[i]] <- -solve(ML$hessian)

      if(progress == TRUE)
        cat("   Time elapsed:", (proc.time()-ptm)[3], "s", "\n")
    }
  }

  mis_coef  <- as.matrix(coef[(p1 + 1):(p1 + p2), ])
  if(p2 == 1){
    mis_coef  <- t(mis_coef)
    rownames(mis_coef) <- names(mis_model$coefficients)
  }

  colnames(mis_coef) <- paste(rho)
  coef <- as.matrix(coef[1:(p1), ])
  if(p1 == 1){
    coef <- t(coef)
    rownames(coef) <- names(out_model$coefficients)
  }
  colnames(coef) <- paste(rho)

  max.info <- lapply(max.info, setNames, nm = paste(rho))
  ml_object<-list(coef = coef, rho = rho, mis_coef = mis_coef, mis_model = mis_model, out_model = out_model, Xz = Xz,
              Xy = Xy, y = y, z = z, value = value, vcov = sigma, max.info = max.info)
  return(ml_object)

}