Merge pull request #141 from PLN-team/sandwich-var

Fixing Sandwich variance - Issue #140

Author: jchiquet

.github/workflows/R-CMD-check.yaml

@@ -19,7 +19,6 @@ jobs:
       matrix:
         config:
           - {os: macOS-latest,   r: 'release'}
-          - {os: macOS-latest,   r: 'oldrel-1'}
           - {os: windows-latest, r: 'release'}
           - {os: windows-latest, r: 'oldrel-1'}
           - {os: ubuntu-latest,  r: 'release'}

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: PLNmodels
 Title: Poisson Lognormal Models
-Version: 1.2.1
+Version: 1.2.2
 Authors@R: c(
     person("Julien", "Chiquet", , "julien.chiquet@inrae.fr", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-3629-3429")),

NEWS.md

@@ -1,5 +1,6 @@
 # Current version
 
+* fix sandwich variance estimation (PR #140)
 * fix use of native pipe to ensure compatibility with R 3.6 (merge PR #125, fix #124)
 
 # PLNmodels 1.2.0 (2024-03-05)

R/PLNfit-S3methods.R

@@ -156,7 +156,7 @@ sigma.PLNfit <- function(object, ...) {
 #' @description Extracts univariate standard errors for the estimated coefficient of B. Standard errors are computed from the (approximate) Fisher information matrix.
 #'
 #' @param object an R6 object with class PLNfit
-#' @param type string describing the type of variance approximation: "variational", "jackknife", "sandwich" (only for fixed covariance). Default is "variational".
+#' @param type string describing the type of variance approximation: "variational", "jackknife", "sandwich". Default is "sandwich".
 #' @param parameter string describing the target parameter: either B (regression coefficients) or Omega (inverse residual covariance)
 #'
 #' @seealso [vcov.PLNfit()] for the complete variance covariance estimation of the coefficient
@@ -166,16 +166,16 @@ sigma.PLNfit <- function(object, ...) {
 #' data(trichoptera)
 #' trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
 #' myPLN <- PLN(Abundance ~ 1 + offset(log(Offset)), data = trichoptera,
-#'               control = PLN_param(config_post = list(variational_var = TRUE)))
+#'               control = PLN_param(config_post = list(sandwich_var = TRUE)))
 #' standard_error(myPLN)
 #' @export
-standard_error <- function(object, type = c("variational", "jackknife", "sandwich"), parameter = c("B", "Omega")) {
+standard_error <- function(object, type = c("sandwich", "variational", "jackknife"), parameter = c("B", "Omega")) {
   UseMethod("standard_error", object)
 }
 
 #' @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("sandwich", "variational", "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("sandwich", "variational", "jackknife", "bootstrap"), 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()
 }

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,33 +233,32 @@ 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) {
+      vcov_sand <- get_sandwich_variance_B(Y, X, private$A,
+                                           private$S, private$Sigma, diag(private$Omega)
+      )
+      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],
                      X = X[-i, , drop = FALSE],
                      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(`+`) %>%
@@ -298,21 +286,17 @@ PLNfit <- R6Class(
                      O = O[resample, , drop = FALSE],
                      w = w[resample])
         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)
+          data   <- lapply(data, torch_tensor, device = config$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 +439,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 +480,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 +909,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,39 +948,14 @@ 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(`+`) %>%
         `dimnames<-`(dimnames(private$B))
       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(

R/RcppExports.R

@@ -89,3 +89,7 @@ cpp_test_packing <- function() {
     .Call('_PLNmodels_cpp_test_packing', PACKAGE = 'PLNmodels')
 }
 
+get_sandwich_variance_B <- function(Y, X, A, S, Sigma, Diag_Omega) {
+    .Call('_PLNmodels_get_sandwich_variance_B', PACKAGE = 'PLNmodels', Y, X, A, S, Sigma, Diag_Omega)
+}
+

R/utils.R

@@ -55,23 +55,26 @@ config_post_default_PLNLDA <-
     jackknife       = FALSE,
     bootstrap       = 0L,
     rsquared        = TRUE,
-    variational_var = FALSE
+    variational_var = FALSE,
+    sandwich_var    = FALSE
   )
 
 config_post_default_PLNPCA <-
   list(
     jackknife       = FALSE,
     bootstrap       = 0L,
     rsquared        = TRUE,
-    variational_var = FALSE
+    variational_var = FALSE,
+    sandwich_var    = FALSE
   )
 
 config_post_default_PLNmixture <-
   list(
     jackknife       = FALSE,
     bootstrap       = 0L,
     rsquared        = TRUE,
-    variational_var = FALSE
+    variational_var = FALSE,
+    sandwich_var    = FALSE
   )
 
 status_to_message <- function(status) {

inst/check/variance_estimates.R

@@ -1,37 +1,68 @@
 library(tidyverse)
 library(PLNmodels)
+library(future.apply)
 set.seed(1234)
 
-nb_cores <- 10
-options(future.fork.enable = TRUE)
+plan(multisession, workers = 10)
 
-params <- PLNmodels:::create_parameters(n = 50, p = 10, d = 1, depths = 1e3)
+rmse <- function(theta_hat, theta_star) {
+  sqrt(mean((theta_hat - theta_star)^2))
+}
+
+params <- PLNmodels:::create_parameters(n = 250, p = 10, d = 1, depths = 1e3)
 X <- params$X
 B <- params$B
-Y <- rPLN(n = nrow(X), mu = X %*% B, Sigma = params$Sigma, depths = params$depths)
+conf <- list(variational_var = TRUE, jackknife = FALSE, bootstrap = 50, sandwich_var = TRUE)
+seq_n <- c(50, 100, 250)
+
+one_simu <- function(s) {
+  cat("\nsimu", s)
+  do.call(rbind, lapply(seq_n, function(n) {
+    X_sub <- X[1:n, , drop = FALSE]
+    Y <- rPLN(n = n, mu = X_sub %*% B, Sigma = params$Sigma, depths = params$depths)
+    data <- prepare_data(Y, X_sub, offset = "none")
+    logO <- attr(Y, "offsets")
+    model <- PLN(Abundance ~ 0 + . + offset(logO), data = data, control = PLN_param(trace = FALSE, config_post = conf))
+
+    B_hat <- coef(model)
+    vcov_sandwich    <- attr(coef(model), "vcov_sandwich")
+    vcov_bootstrap   <- attr(coef(model), "vcov_bootstrap")
+    vcov_variational <- attr(coef(model), "vcov_variational")
+
+    data.frame(rmse = rmse(B_hat, B),
+               cover_sandwich    = mean(abs(as.numeric(B_hat - B) %*% solve(chol(vcov_sandwich))) < 1.96),
+               cover_bootstrap   = mean(abs(as.numeric(B_hat - B) %*% solve(chol(vcov_bootstrap))) < 1.96),
+               cover_variational = mean(abs(as.numeric(B_hat - B) %*% solve(chol(vcov_variational))) < 1.96),
+               sample_size = n,
+               simu = s)
+  }))
+}
 
+res <- do.call(rbind, lapply(1:50, one_simu))
+
+res %>%
+  pivot_longer(-c(rmse, sample_size, simu), values_to = "coverage", names_to = "estimator") %>%
+  ggplot() + aes(x = estimator, y = coverage, fill = factor(sample_size)) + geom_boxplot()
+
+### Single test
+Y <- rPLN(n = nrow(X), mu = X %*% B, Sigma = params$Sigma, depths = params$depths)
 data <- prepare_data(Y, X, offset = "none")
 logO <- attr(Y, "offsets")
-
-conf <- list(variational_var = TRUE, jackknife = TRUE, bootstrap = nrow(Y))
-future::plan("multicore", workers = nb_cores)
 model <- PLN(Abundance ~ 0 + . + offset(logO), data = data, control = PLN_param(config_post = conf))
-future::plan("sequential")
 
 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_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))) ) %>%
+  se = c(B_se_var, B_se_bt, B_se_sw),
+  method = rep(c("variational", "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) +
   geom_abline(slope = 1, intercept = 0) + labs(x = "True value", y = "Mean estimate") + theme_bw() -> p
 
 print(p)
-