From 4671672cf2b72127ce38cc25f51686e575643764 Mon Sep 17 00:00:00 2001 From: Walmes Zeviani <walmes@ufpr.br> Date: Thu, 7 Jan 2016 15:33:17 -0200 Subject: [PATCH] Add tags and tidy code and text - Add @title, @name and @author. - Tidy code and text. --- R/mc_bias_correct_std.R | 65 ++++++++++++++++++++++------------------- 1 file changed, 35 insertions(+), 30 deletions(-) diff --git a/R/mc_bias_correct_std.R b/R/mc_bias_correct_std.R index c9c5c16..5cd1b4a 100644 --- a/R/mc_bias_correct_std.R +++ b/R/mc_bias_correct_std.R @@ -1,44 +1,49 @@ -#' Bias-corrected standard error for regression parameters +#' @title Bias-corrected standard error for regression parameters +#' @name mc_bias_corrected_std +#' @author Wagner Hugo Bonat #' -#' @description Compute bias-corrected standard error for regression parameters in the context -#' of clustered observations. It is also robust and has improved finite sample properties. +#' @description Compute bias-corrected standard error for regression +#' parameters in the context of clustered observations. It is also +#' robust and has improved finite sample properties. #' #' @param object An object of mcglm class. -#' @param id a vector which identifies the clusters. The length and order of id should be the -#' same as the number of observations. Data are assumed to be sorted so that observations on a cluster -#' are contiguous rows for all entities in the formula. -#' @return A matrix. Note that the function assumes that the data are in the correct order. +#' @param id a vector which identifies the clusters. The length and +#' order of id should be the same as the number of +#' observations. Data are assumed to be sorted so that observations +#' on a cluster are contiguous rows for all entities in the formula. +#' @return A matrix. Note that the function assumes that the data are in +#' the correct order. #' @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 <- 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]) })) - 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]] + 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(temp_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)) + 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