From 59e71b95efa4c50a4bbe1de64639d7a30bb1753e Mon Sep 17 00:00:00 2001
From: wbonat <wbonat@gmail.com>
Date: Tue, 6 Oct 2015 15:32:21 +0200
Subject: [PATCH] Fix matrix H
---
R/mc_bias_correct_std.R | 37 +++++++++++++++++++++++++++----------
1 file changed, 27 insertions(+), 10 deletions(-)
diff --git a/R/mc_bias_correct_std.R b/R/mc_bias_correct_std.R
index cfe7f9f..718e9c1 100644
--- a/R/mc_bias_correct_std.R
+++ b/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)
}
+
--
GitLab