Fix matrix H

R/mc_bias_correct_std.R

@@ -11,17 +11,34 @@
 #' @export
 
 mc_bias_corrected_std <- function(object, id) {
-    inv_M <- object$inv_S_beta
-    temp_data <- data.frame(res = object$residuals, id)
-    temp_data_group <- split(temp_data, temp_data$id)
-    D <- bdiag(lapply(object$mu_list, function(x) x$D))
-    r_rT <- bdiag(lapply(temp_data_group, function(x) {
+  inv_M <- -object$inv_S_beta
+  temp_data <- data.frame(res = object$residuals, id)
+  temp_data_group <- split(temp_data, temp_data$id)
+  D <- bdiag(lapply(object$mu_list, function(x) x$D))
+  R <- bdiag(lapply(temp_data_group, function(x) {
         tcrossprod(x[, 1])
     }))
-    tD_invC <- t(D) %*% object$inv_C
-    H <- Matrix(D %*% inv_M %*% tD_invC, sparse = TRUE)
-    IH <- Diagonal(object$n_obs, 1) - H
-    inv_IH <- solve(IH)
-    output <- sqrt(diag(inv_M %*% tD_invC %*% inv_IH %*% r_rT %*% inv_IH %*% t(tD_invC) %*% inv_M))
+  uni_id <- unique(id)
+  n_id <- length(uni_id)
+  Hi <- list()
+  Di <- list()
+  inv_Ci <- list()
+  for (i in 1:n_id) {
+    idTF <- id == uni_id[i]
+    if(sum(idTF) == 1) {
+      Di[[i]] <- Matrix(D[idTF,], nrow = sum(idTF), ncol = dim(D)[2])
+      } else {
+        Di[[i]] <- D[idTF,]
+      }
+    inv_Ci[[i]] <- object$inv_C[idTF,idTF]
+    Hi[[i]] <- Di[[i]]%*%inv_M%*%t(Di[[i]])%*%inv_Ci[[i]]
+    }
+    H <- bdiag(Hi)
+    I <- Diagonal(dim(data)[1], 1)
+    inv_IH <- solve(I - H)
+    Vbeta = inv_M%*%(t(D)%*%object$inv_C%*%inv_IH%*%R%*%inv_IH%*%
+                       object$inv_C%*%D)%*%inv_M
+    output = sqrt(diag(Vbeta))
     return(output)
 }
+