R/louisvcov.R

In glmfitmiss: Fitting GLMs with Missing Data in Both Responses and Covariates

#---------Computing the variance using Louis(1982) ---------------------
# ------input data should be having a weight variable wgt -------------------

#' @importFrom stats model.matrix vcov
#' @importFrom MASS ginv
#' @importFrom utils globalVariables
#' @importFrom data.table as.data.table := .I setorder
#' @keywords internal
louisvcov <- function (formula, data=final_df, family=binomial, correct="brglmFit"){

  #predictors of the formula, the scope of the variable resp is defined as global
  predictor_vars <- attr(terms(formula), "term.labels")
  p <- length(predictor_vars)+1
  resp <- all.vars(formula)[1]
  predictor_data <- model.matrix(terms(formula), data)
  if(correct=="brglmFit"){
    fit_y<-suppressWarnings(stats::glm(formula, family = family, data = data, weights=wgt, method = brglm2::brglmFit))
    vcov_beta <- stats::vcov(fit_y)
    pi<-stats::fitted(fit_y)
    p1=as.vector(pi*(1-pi))
    wvec<-diag(p1)
    X <- model.matrix(fit_y)
    W_sqrt <- sqrt(diag(p1))
    H <- W_sqrt %*% X %*% vcov_beta %*% t(X) %*% W_sqrt
    hi <- diag(H)
    # hi <- cooks.distance(fit_y) #alternate of computing hi from H
    # H computation takes longer time for larger data, cooks distance is recommended
    resp_score_fn <- data[, resp] - pi+ hi*(0.5-pi)
  }
  else{
    fit_y<-suppressWarnings(stats::glm(formula, family = family, data = data, weights=wgt))
    vcov_beta <- stats::vcov(fit_y)
    pi<-stats::fitted(fit_y)
    resp_score_fn <- data[, resp] - pi
  }

  wij <- as.numeric(data$wgt)
  # ------------first part: negQdd (second derivative of Q)-------------
  negQdd<-MASS::ginv(vcov_beta)

  # -----second part-------------------------------------------------------------
  part2 <- matrix(0, ncol = p, nrow = p)
  data_len <- nrow(data)

  for (i in 1:data_len) {
    t1 <- predictor_data[i,]* (resp_score_fn[i] * wij[i])
    t2 <- predictor_data[i,]* resp_score_fn[i]
    tmp <-  t1%*% t(t2)
    part2 <- part2 + tmp
  }

  # -----third part-------------------------------------------------------------

  weightsTemp <- resp_score_fn * wij
  t1 <- as.data.frame(sweep(predictor_data,1, weightsTemp, '*'))
  # t1 <- as.data.frame(predictor_data * (resp_score_fn * wij))

  groups <- unique(data$grp)

  # Create an empty data frame to store the column sums
  sumsbygrp <- data.frame(grp = groups)

  # Calculate the column sums for each group
  for (var_name in names(t1)) {
    sumsbygrp[var_name] <- tapply(t1[[var_name]], data$grp, sum)
  }

  #removing the grp variable to avoid error in matrix multiplication
  sumsbygrp$grp <- NULL

  part3 <- matrix(0, ncol = p, nrow = p)

  for (i in 1:length(groups)) {
    q11 <- as.matrix(sumsbygrp[i, ])
    q12 <- t(q11) %*% q11
    part3 <- part3 + q12
  }

  # -----final part-------------------------------------------------------------
  Info<-negQdd - (part2-part3)
  invInfo<-MASS::ginv(Info)
  vcov_beta<-invInfo[1:p, 1:p]

  return(vcov_beta)

}