BIOST561FinalProject

Zihang Yu

2024-06-07

Warning: this is a demo package for vcovCR.glmerMod, the complete version is on the way.

Introduction

The goal is to write a function which will provide robust variances for generalized linear mixed models fit using the glmer function in R.

vcovCR.glmerModZY <- function(obj, cluster, type = "classic") {
  # put some check here to ensure the input is correct
  # If cluster not specified, will be set to attr(obj,"cluster")
  # if (is.null(cluster)) {
  #   cluster <- attr(obj, "cluster")
  # }

  link <- family(obj)$link

  beta <- matrix(fixef(obj), ncol = 1)

  np <- length(beta)

  gamma <- matrix(unlist(ranef(obj)), ncol = 1)

  nq <- length(gamma)

  X <- model.matrix(obj,type = "fixed")
  Z <- model.matrix(obj, type = "random")  # Z matrix for random effects
  Y <- obj@resp$y

  eta <- predict(obj, type = "link")
  ginv_eta <- predict(obj, type = "response")

  if (link == "identity"){
    delta = diag(nobs(obj)) # diag matrix
    deltainv = diag(1/diag(delta))# a more efficient way to get solve(delta)
  }

  else if (link == "logit"){
    delta <- diag(exp(eta)/(1+exp(eta))^2)
    deltainv <- diag(1/diag(delta))
  }

  else if (link == "log") {
    delta <- diag(exp(eta))
    deltainv <- diag(1/diag(delta))
  }

  P <- deltainv %*% (Y - ginv_eta) + eta
  e <- matrix(P - X %*% beta, ncol = 1)
  XtVX <- vcov(obj)  # model based variance
  #theta <- as.data.frame(VarCorr(obj)) # This is where the problem is, when you fit different models, this is different

  sigma2 <- sigma(obj)^2
  lambda <- getME(obj,"Lambda")
  R <- lambda %*% t(lambda) * sigma2

  G <- ngrps(obj)["fcluster"]
  sum <- matrix(0, np, np)

  for (g in 1:G){
    grp = ctdata$fcluster == g
    ng = sum(grp)

    # V = Z[grp,]%*%R%*%t(Z[grp,]) + deltainv[grp,grp]%*%Sigma%*%deltainv[grp,grp]
    # Vinv = solve(V)
    Sigma <- sigma2 * diag(ng)
    WB_A <- diag(1/diag(deltainv[grp,grp]%*%Sigma%*%deltainv[grp,grp]))

    WB_U <- Z[grp,]
    WB_C <- solve(R)
    WB_V <- t(Z[grp,])

    W <- WoodburyMatrix(A = WB_A, U = WB_U, B = WB_C, V = WB_V)
    Vinv <- solve(W)

    H = X[grp, ] %*% XtVX %*% t(X[grp, ]) %*% Vinv
    Q = t(X[grp, ]) %*% Vinv %*% X[grp, ] %*% XtVX

    # if loop, choose A, F and c based on robust variance form specified by the input 'type'
    F = diag(ng)
    A = diag(np)

    sum = sum + A %*% t(X[grp, ]) %*% Vinv %*% t(F) %*% e[grp, ] %*% t(e[grp, ]) %*% F %*% Vinv %*% X[grp, ] %*% A
  }

  c = 1
  robustVar <- c * XtVX %*% sum %*% XtVX
  diag(robustVar)  # Return diagonal elements representing the variances
}

Example

library("sandwich")
library("clubSandwich")

## Registered S3 method overwritten by 'clubSandwich':
##   method    from    
##   bread.mlm sandwich

library("lme4")

## Loading required package: Matrix

library("WoodburyMatrix")

# # Import data and fit lmer model, creat factors for clustering
# url <- "https://github.com/ZAYNYU/BIOST561FinalProject/blob/main/vignettes/ctdata.Rdata"  # Use the download link
# temp <- tempfile()  # Creates a temporary file which will be deleted later
# download.file(url, temp, mode = "wb")  # Download the file in binary mode
# load(temp)  # Load the .RData file into your R session
# unlink(temp)  # Remove the temporary file immediately after loading
# load(url("https://github.com/ZAYNYU/BIOST561FinalProject/blob/main/vignettes/ctdata.Rdata"))
load(url("https://raw.githubusercontent.com/ZAYNYU/BIOST561FinalProject/main/vignettes/ctdata.Rdata"))

ctdata$fcluster = factor(ctdata$hdist)
ctdata$ftime = factor(ctdata$time)
ctdata$clustime = interaction(ctdata$fcluster,ctdata$ftime)
rslt = lmer(ct ~ ftime + treat + (1 | fcluster) + (1 | clustime), data=ctdata)
vcovCR.glmerModZY(rslt, ctdata$fcluster)

##  (Intercept)       ftime1       ftime2       ftime3       ftime4        treat 
## 0.0001379472 0.0002188822 0.0002572821 0.0003607838 0.0003631747 0.0001417230