From 2e915c6e8bab9b31036ba9045cf3302cdef9d3dc Mon Sep 17 00:00:00 2001 From: Julien Chiquet Date: Wed, 20 Nov 2024 18:13:29 +0100 Subject: [PATCH 1/7] putting sandwich back to business --- R/PLNfit-S3methods.R | 10 ++- R/PLNfit-class.R | 108 +++++++++------------------ R/utils-zipln.R | 2 +- R/utils.R | 9 ++- inst/check/variance_estimates.R | 13 ++-- man/PLNfit.Rd | 2 +- man/PLNfit_fixedcov.Rd | 50 +------------ man/PLNnetworkfit.Rd | 2 +- man/standard_error.Rd | 4 +- tests/testthat/test-standard-error.R | 13 +++- 10 files changed, 70 insertions(+), 143 deletions(-) diff --git a/R/PLNfit-S3methods.R b/R/PLNfit-S3methods.R index 8fd95053..7dfe0dfd 100644 --- a/R/PLNfit-S3methods.R +++ b/R/PLNfit-S3methods.R @@ -175,7 +175,7 @@ standard_error <- function(object, type = c("variational", "jackknife", "sandwic #' @describeIn standard_error Component-wise standard errors of B in [`PLNfit`] #' @export -standard_error.PLNfit <- function(object, type = c("variational", "jackknife", "bootstrap", "sandwich"), parameter = c("B", "Omega")) { +standard_error.PLNfit <- function(object, type = c("variational", "sandwich", "jackknife", "bootstrap"), parameter = c("B", "Omega")) { type <- match.arg(type) par <- match.arg(parameter) if (type == "variational" & is.null(attr(object$model_par$B, "variance_variational"))) @@ -184,14 +184,14 @@ standard_error.PLNfit <- function(object, type = c("variational", "jackknife", " stop("Jackknife estimation not available: rerun by setting `jackknife = TRUE` in the control list.") if (type == "bootstrap" & is.null(attr(object$model_par$B, "variance_bootstrap"))) stop("Bootstrap estimation not available: rerun by setting `bootstrap > 0` in the control list.") - if (type == "sandwich") - stop("Sandwich estimator is only available for fixed covariance / precision matrix.") + if (type == "sandwich" & is.null(attr(object$model_par$B, "variance_sandwich"))) + stop("Sandwich estimator not available: rerun by setting `sandwich_var = TRUE` in the control list.") attr(object$model_par[[par]], paste0("variance_", type)) %>% sqrt() } #' @describeIn standard_error Component-wise standard errors of B in [`PLNfit_fixedcov`] #' @export -standard_error.PLNfit_fixedcov <- function(object, type = c("variational", "jackknife", "bootstrap", "sandwich"), parameter = c("B", "Omega")) { +standard_error.PLNfit_fixedcov <- function(object, type = c("variational", "sandwich", "jackknife", "bootstrap", "sandwich"), parameter = c("B", "Omega")) { type <- match.arg(type) par <- match.arg(parameter) if (par == "Omega") @@ -202,5 +202,7 @@ standard_error.PLNfit_fixedcov <- function(object, type = c("variational", "jack stop("Jackknife estimation not available: rerun by setting `jackknife = TRUE` in the control list.") if (type == "bootstrap" & is.null(attr(object$model_par$B, "variance_bootstrap"))) stop("Bootstrap estimation not available: rerun by setting `bootstrap > 0` in the control list.") + if (type == "sandwich" & is.null(attr(object$model_par$B, "variance_sandwich"))) + stop("Sandwich estimator not available: rerun by setting `sandwich_var = TRUE` in the control list.") attr(object$model_par[[par]], paste0("variance_", type)) %>% sqrt() } diff --git a/R/PLNfit-class.R b/R/PLNfit-class.R index b964cab5..66e61b85 100644 --- a/R/PLNfit-class.R +++ b/R/PLNfit-class.R @@ -191,7 +191,7 @@ PLNfit <- R6Class( variance_variational = function(X, config = config_default_nlopt) { ## Variance of B for n data points fisher <- Matrix::bdiag(lapply(1:self$p, function(j) { - crossprod(X, private$A[, j] * X) # t(X) %*% diag(A[, i]) %*% X + crossprod(X, private$A[, j] * X) # t(X) %*% diag(A[, j]) %*% X })) vcov_B <- tryCatch(Matrix::solve(fisher), error = function(e) {e}) if (is(vcov_B, "error")) { @@ -220,22 +220,11 @@ PLNfit <- R6Class( compute_vcov_from_resamples = function(resamples){ B_list = resamples %>% map("B") - #print (B_list) vcov_B = lapply(seq(1, ncol(private$B)), function(B_col){ param_ests_for_col = B_list %>% map(~.x[, B_col]) param_ests_for_col = do.call(rbind, param_ests_for_col) - #print (param_ests_for_col) row_vcov = cov(param_ests_for_col) }) - #print ("vcov blocks") - #print (vcov_B) - - #B_vcov <- resamples %>% map("B") %>% map(~( . )) %>% reduce(cov) - - #var_jack <- jacks %>% map("B") %>% map(~( (. - B_jack)^2)) %>% reduce(`+`) %>% - # `dimnames<-`(dimnames(private$B)) - #B_hat <- private$B[,] ## strips attributes while preserving names - vcov_B = Matrix::bdiag(vcov_B) %>% as.matrix() rownames(vcov_B) <- colnames(vcov_B) <- @@ -244,18 +233,35 @@ PLNfit <- R6Class( ## Hack to make sure that species is first and varies slowest apply(1, paste0, collapse = "_") - #print (pheatmap::pheatmap(vcov_B, cluster_rows=FALSE, cluster_cols=FALSE)) - - - #names = lapply(bootstrapped_df$cov_mat, function(m){ colnames(m)}) %>% unlist() - #rownames(bootstrapped_vhat) = names - #colnames(bootstrapped_vhat) = names - vcov_B = methods::as(vcov_B, "dgCMatrix") return(vcov_B) }, + vcov_sandwich_B = function(Y, X) { + getMat_iCnB <- function(i) { + a_i <- as.numeric(private$A[i, ]) + s2_i <- as.numeric(private$S[i, ]**2) + omega <- as.numeric(diag(private$Omega)) + diag_mat_i <- diag(1/a_i + s2_i^2 / (1 + s2_i * (a_i + omega))) + solve(private$Sigma + diag_mat_i) + } + YmA <- Y - private$A + Dn <- matrix(0, self$d*self$p, self$d*self$p) + Cn <- matrix(0, self$d*self$p, self$d*self$p) + for (i in 1:self$n) { + xxt_i <- tcrossprod(X[i, ]) + Cn <- Cn - kronecker(getMat_iCnB(i) , xxt_i) / self$n + Dn <- Dn + kronecker(tcrossprod(YmA[i,]), xxt_i) / self$n + } + Cn_inv <- solve(Cn) + dim_names <- dimnames(attr(private$B, "vcov_variational")) + vcov_sand <- ((Cn_inv %*% Dn %*% Cn_inv) / self$n) %>% `dimnames<-`(dim_names) + attr(private$B, "vcov_sandwich") <- vcov_sand + attr(private$B, "variance_sandwich") <- matrix(diag(vcov_sand), nrow = self$d, ncol = self$p, + dimnames = dimnames(private$B)) + }, + variance_jackknife = function(Y, X, O, w, config = config_default_nlopt) { jacks <- future.apply::future_lapply(seq_len(self$n), function(i) { data <- list(Y = Y[-i, , drop = FALSE], @@ -263,14 +269,11 @@ PLNfit <- R6Class( O = O[-i, , drop = FALSE], w = w[-i]) args <- list(data = data, - # params = list(B = private$B, - # M = matrix(0, self$n-1, self$p), - # S = private$S[-i, , drop = FALSE]), params = do.call(compute_PLN_starting_point, data), config = config) optim_out <- do.call(private$optimizer$main, args) optim_out[c("B", "Omega")] - }) + }, future.seed = TRUE, future.scheduling = structure(TRUE, ordering = "random")) B_jack <- jacks %>% map("B") %>% reduce(`+`) / self$n var_jack <- jacks %>% map("B") %>% map(~( (. - B_jack)^2)) %>% reduce(`+`) %>% @@ -300,19 +303,15 @@ PLNfit <- R6Class( if (config$backend == "torch") # Convert data to torch tensors data <- lapply(data, torch_tensor, device = config$device) # list with Y, X, O, w - #print (data$Y$device) - args <- list(data = data, - # params = list(B = private$B, M = matrix(0,self$n,self$p), S = private$S[resample, ]), params = do.call(compute_PLN_starting_point, data), config = config) if (config$backend == "torch") # Convert data to torch tensors args$params <- lapply(args$params, torch_tensor, requires_grad = TRUE, device = config$device) # list with B, M, S optim_out <- do.call(private$optimizer$main, args) - #print (optim_out) optim_out[c("B", "Omega", "monitoring")] - }) + }, future.seed = TRUE, future.scheduling = structure(TRUE, ordering = "random")) B_boots <- boots %>% map("B") %>% reduce(`+`) / n_resamples attr(private$B, "variance_bootstrap") <- @@ -455,7 +454,7 @@ PLNfit <- R6Class( #' * jackknife boolean indicating whether jackknife should be performed to evaluate bias and variance of the model parameters. Default is FALSE. #' * bootstrap integer indicating the number of bootstrap resamples generated to evaluate the variance of the model parameters. Default is 0 (inactivated). #' * variational_var boolean indicating whether variational Fisher information matrix should be computed to estimate the variance of the model parameters (highly underestimated). Default is FALSE. - #' * rsquared boolean indicating whether approximation of R2 based on deviance should be computed. Default is TRUE + #' * sandwich_var boolean indicating whether sandwich estimator should be computed to estimate the variance of the model parameters (highly underestimated). Default is FALSE. #' * trace integer for verbosity. should be > 1 to see output in post-treatments postTreatment = function(responses, covariates, offsets, weights = rep(1, nrow(responses)), config_post, config_optim, nullModel = NULL) { ## PARAMATERS DIMNAMES @@ -496,6 +495,11 @@ PLNfit <- R6Class( } private$variance_bootstrap(responses, covariates, offsets, weights, n_resamples=config_post$bootstrap, config = config_optim) } + ## 5. compute and store matrix of standard variances for B with sandwich correction approximation + if (config_post$sandwich_var) { + if(config_post$trace > 1) cat("\n\tComputing sandwich estimator of the variance...") + private$vcov_sandwich_B(responses, covariates) + } }, #' @description Predict position, scores or observations of new data. @@ -920,25 +924,8 @@ PLNfit_fixedcov <- R6Class( optim_out <- do.call(private$optimizer$main, args) do.call(self$update, optim_out) private$Sigma <- solve(optim_out$Omega) - }, - - #' @description Update R2, fisher and std_err fields after optimization - #' @param config_post a list for controlling the post-treatments (optional bootstrap, jackknife, R2, etc.). See details - #' @param config_optim a list for controlling the optimization parameter. See details - #' @details The list of parameters `config` controls the post-treatment processing, with the following entries: - #' * trace integer for verbosity. should be > 1 to see output in post-treatments - #' * jackknife boolean indicating whether jackknife should be performed to evaluate bias and variance of the model parameters. Default is FALSE. - #' * bootstrap integer indicating the number of bootstrap resamples generated to evaluate the variance of the model parameters. Default is 0 (inactivated). - #' * variational_var boolean indicating whether variational Fisher information matrix should be computed to estimate the variance of the model parameters (highly underestimated). Default is FALSE. - #' * rsquared boolean indicating whether approximation of R2 based on deviance should be computed. Default is TRUE - postTreatment = function(responses, covariates, offsets, weights = rep(1, nrow(responses)), config_post, config_optim, nullModel = NULL) { - super$postTreatment(responses, covariates, offsets, weights, config_post, config_optim, nullModel) - ## 6. compute and store matrix of standard variances for B with sandwich correction approximation - if (config_post$sandwich_var) { - if(config_post$trace > 1) cat("\n\tComputing sandwich estimator of the variance...") - private$vcov_sandwich_B(responses, covariates) - } } + ), private = list( ## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -976,7 +963,7 @@ PLNfit_fixedcov <- R6Class( config = config) optim_out <- do.call(private$optimizer$main, args) optim_out[c("B", "Omega")] - }, future.seed = TRUE) + }, future.seed = TRUE, future.scheduling = structure(TRUE, ordering = "random")) B_jack <- jacks %>% map("B") %>% reduce(`+`) / self$n var_jack <- jacks %>% map("B") %>% map(~( (. - B_jack)^2)) %>% reduce(`+`) %>% @@ -984,31 +971,6 @@ PLNfit_fixedcov <- R6Class( B_hat <- private$B[,] ## strips attributes while preserving names attr(private$B, "bias") <- (self$n - 1) * (B_jack - B_hat) attr(private$B, "variance_jackknife") <- (self$n - 1) / self$n * var_jack - }, - - vcov_sandwich_B = function(Y, X) { - getMat_iCnB <- function(i) { - a_i <- as.numeric(private$A[i, ]) - s2_i <- as.numeric(private$S[i, ]**2) - # omega <- as.numeric(1/diag(private$Sigma)) - # diag_mat_i <- diag(1/a_i + s2_i^2 / (1 + s2_i * (a_i + omega))) - diag_mat_i <- diag(1/a_i + .5 * s2_i^2) - solve(private$Sigma + diag_mat_i) - } - YmA <- Y - private$A - Dn <- matrix(0, self$d*self$p, self$d*self$p) - Cn <- matrix(0, self$d*self$p, self$d*self$p) - for (i in 1:self$n) { - xxt_i <- tcrossprod(X[i, ]) - Cn <- Cn - kronecker(getMat_iCnB(i) , xxt_i) / (self$n) - Dn <- Dn + kronecker(tcrossprod(YmA[i,]), xxt_i) / (self$n) - } - Cn_inv <- solve(Cn) - dim_names <- dimnames(attr(private$B, "vcov_variational")) - vcov_sand <- ((Cn_inv %*% Dn %*% Cn_inv) / self$n) %>% `dimnames<-`(dim_names) - attr(private$B, "vcov_sandwich") <- vcov_sand - attr(private$B, "variance_sandwich") <- matrix(diag(vcov_sand), nrow = self$d, ncol = self$p, - dimnames = dimnames(private$B)) } ), active = list( diff --git a/R/utils-zipln.R b/R/utils-zipln.R index 59e9233a..f62cdc0f 100644 --- a/R/utils-zipln.R +++ b/R/utils-zipln.R @@ -27,7 +27,7 @@ extract_model_zi <- function(call, envir) { call_args <- c(as.list(call_args), list(xlev = attr(call$formula, "xlevels"), na.action = NULL)) ## Extract terms for ZI and PLN components - terms <- .extract_terms_zi(as.formula(eval(call$formula, env = envir))) + terms <- .extract_terms_zi(as.formula(eval(call$formula, envir = envir))) ## eval the call in the parent environment with adjustement due to ZI terms call_args$formula <- terms$formula frame <- do.call(stats::model.frame, call_args, envir = envir) diff --git a/R/utils.R b/R/utils.R index 2ccba19c..204861b8 100644 --- a/R/utils.R +++ b/R/utils.R @@ -55,7 +55,8 @@ config_post_default_PLNLDA <- jackknife = FALSE, bootstrap = 0L, rsquared = TRUE, - variational_var = FALSE + variational_var = FALSE, + sandwich_var = FALSE ) config_post_default_PLNPCA <- @@ -63,7 +64,8 @@ config_post_default_PLNPCA <- jackknife = FALSE, bootstrap = 0L, rsquared = TRUE, - variational_var = FALSE + variational_var = FALSE, + sandwich_var = FALSE ) config_post_default_PLNmixture <- @@ -71,7 +73,8 @@ config_post_default_PLNmixture <- jackknife = FALSE, bootstrap = 0L, rsquared = TRUE, - variational_var = FALSE + variational_var = FALSE, + sandwich_var = FALSE ) status_to_message <- function(status) { diff --git a/inst/check/variance_estimates.R b/inst/check/variance_estimates.R index 65cafab7..8eebc43c 100644 --- a/inst/check/variance_estimates.R +++ b/inst/check/variance_estimates.R @@ -13,7 +13,7 @@ Y <- rPLN(n = nrow(X), mu = X %*% B, Sigma = params$Sigma, depths = params$depth data <- prepare_data(Y, X, offset = "none") logO <- attr(Y, "offsets") -conf <- list(variational_var = TRUE, jackknife = TRUE, bootstrap = nrow(Y)) +conf <- list(variational_var = TRUE, jackknife = TRUE, bootstrap = nrow(Y), sandwich_var = TRUE) future::plan("multicore", workers = nb_cores) model <- PLN(Abundance ~ 0 + . + offset(logO), data = data, control = PLN_param(config_post = conf)) future::plan("sequential") @@ -22,12 +22,15 @@ B_hat <- coef(model) B_se_var <- standard_error(model, "variational") B_se_jk <- standard_error(model, "jackknife") B_se_bt <- standard_error(model, "bootstrap") +B_se_var <- standard_error(model, "variational") +B_se_jk <- standard_error(model, "jackknife") +B_se_sw <- standard_error(model, "sandwich") data.frame( - B = rep(c(B), 3), - B_hat = rep(c(B_hat), 3), - se = c(B_se_var, B_se_jk, B_se_bt), - method = rep(c("variational", "jackknife", "bootstrap"), each = length(c(B))) ) %>% + B = rep(c(B), 4), + B_hat = rep(c(B_hat), 4), + se = c(B_se_var, B_se_jk, B_se_bt, B_se_sw), + method = rep(c("variational", "jackknife", "bootstrap", "sandwich"), each = length(c(B))) ) %>% ggplot(aes(x = B, y = B_hat)) + geom_errorbar(aes(ymin = B_hat - 2 * se, ymax = B_hat + 2 * se), color = "blue") + facet_wrap(~ method) + diff --git a/man/PLNfit.Rd b/man/PLNfit.Rd index 303c21cb..f99471c0 100644 --- a/man/PLNfit.Rd +++ b/man/PLNfit.Rd @@ -273,7 +273,7 @@ The list of parameters \code{config} controls the post-treatment processing, wit \item jackknife boolean indicating whether jackknife should be performed to evaluate bias and variance of the model parameters. Default is FALSE. \item bootstrap integer indicating the number of bootstrap resamples generated to evaluate the variance of the model parameters. Default is 0 (inactivated). \item variational_var boolean indicating whether variational Fisher information matrix should be computed to estimate the variance of the model parameters (highly underestimated). Default is FALSE. -\item rsquared boolean indicating whether approximation of R2 based on deviance should be computed. Default is TRUE +\item sandwich_var boolean indicating whether sandwich estimator should be computed to estimate the variance of the model parameters (highly underestimated). Default is FALSE. \item trace integer for verbosity. should be > 1 to see output in post-treatments } } diff --git a/man/PLNfit_fixedcov.Rd b/man/PLNfit_fixedcov.Rd index 52e75903..d0e39325 100644 --- a/man/PLNfit_fixedcov.Rd +++ b/man/PLNfit_fixedcov.Rd @@ -36,7 +36,6 @@ print(myPLN) \itemize{ \item \href{#method-PLNfit_fixedcov-new}{\code{PLNfit_fixedcov$new()}} \item \href{#method-PLNfit_fixedcov-optimize}{\code{PLNfit_fixedcov$optimize()}} -\item \href{#method-PLNfit_fixedcov-postTreatment}{\code{PLNfit_fixedcov$postTreatment()}} \item \href{#method-PLNfit_fixedcov-clone}{\code{PLNfit_fixedcov$clone()}} } } @@ -44,6 +43,7 @@ print(myPLN)
Inherited methods